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PREFACE 

The  contents  of  this  report  consist  of  a  dissertation  written  by  Frederick  M.  Luther 
submitted  in  partial  satisfaction  of  the  requirements  for  the  degree  of  Doctor  of  Philosophy 
in  Engineering-Applied  Science  in  the  Graduate  Division  of  the  University  of  California, 
Davis,  California. 

Dr.  L.  0.  Myrup  of  the  University  of  California  at  Davis  served  as  Luther's  dissertation 
advisor  and  his  interest  in  micro-mesoscale  meteorological  interactions  and  participation  in 
Army  Grant  DA-AMC-28-043-68-G10  are  reflected  in  the  guidance  of  the  work  presented  in  this 
report . 

The  Davis  micrometeorological  data  used  in  testing  the  numerical  model  were  gathered 
during  the  Cooperative  Field  Experiment  in  May  1967  and  processed  under  the  direction  of 
W.  0.  Pruitt,  D.  L.  Morgan,  and  F.  J.  Lourence  of  the  Hydrometeorological  Research  Group  of 
the  Department  of  Water  Science  and  Engineering  at  UC  Davis,  with  the  primary  sponsorship 
coming  from  the  United  States  Army  Electronics  Command,  Atmospheric  Sciences  Laboratory  at 
Fort  Huachuca,  Arizona,  Contract  DA-02-086-AMC-0447(E) .   The  corresponding  upper  air  data 
used  by  Luther  were  gathered  by  members  of  the  ESSA  Weather  Bureau  Research  Station,  Fort 
Huachuca,  Arizona  and  the  Atmospheric  Sciences  Office,  ASL,  USAECOM,  White  Sands  Missile 
Range,  New  Mexico.   Additional  support  for  instrumentation  and  personnel  was  furnished  by 
the  Water  Resources  Center,  University  of  California. 
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ABSTRACT 

A  one-dimensional  numerical  model  is  used  to  compute  temperature  and  moisture  profiles 
and  to  analyze  constituents  of  the  energy  budget  in  the  lower  atmosphere.   The  energy  trans- 
fer processes  involved  in  the  calculation  include  convection,  advection,  absorption  of  solar 
radiation,  exchange  of  terrestrial  radiation,  evaporation,  condensation,  and  conduction  of 
heat  into  the  soil.   The  wind  field  is  specified  for  the  duration  of  the  calculation,  and 
the  mean  vertical  motion  is  assumed  equal  to  zero. 

The  region  of  primary  interest  is  the  first  2000  meters  of  the  atmosphere,  which  is 
the  region  where  the  most  rapid  changes  in  temperature  and  moisture  content  occur.   The 
model  covers  a  region  extending  to  a  height  of  8000  meters  above  the  earth's  surface  in 
order  to  include  the  water  vapor  in  this  region  in  the  terrestrial  radiation  calculation. 
The  model  also  includes  a  layer  of  soil  extending  to  a  depth  of  one  meter,  which  is  needed 
for  the  energy  balance  calculation  at  the  earth's  surface. 

An  expanding  vertical  grid  is  used  in  the  atmospheric  region  in  which  each  grid  length 
is  1.25  times  the  size  of  the  one  below  it,  thus  providing  greater  resolution  in  the  lower 
layers.   The  minimum  grid  length  is  2.0  meters. 

Temperature  and  moisture  fields  are  computed  for  two  twenty-four  hour  periods;  one  be- 
ginning at  1230  LST  on  August  24,  1953  at  O'Neill,  Nebraska  and  the  other  beginning  at 
0100  LST  on  May  4,  1967  at  Davis,  California.   In  both  cases  the  computed  temperature  and 
moisture  fields  agree  well  with  the  observed  values. 

The  O'Neill  calculation  demonstrates  the  ability  of  the  numerical  model  to  predict 
temperature  and  moisture  fields  under  conditions  of  strong  winds.   The  computed  temperature 
profiles  are  shown  to  be  more  accurate  than  those  computed  by  Estoque  for  the  same  twenty- 
four  hour  period. 

The  Davis  calculation  demonstrates  the  ability  of  the  numerical  model  to  predict  tem- 
perature and  moisture  profiles  under  conditions  of  light  winds.   Results  show  that  the  cal- 
culation of  accurate  convective  heat  fluxes  near  the  earth's  surface  is  more  difficult 
under  conditions  of  strong  thermal  stability.   Accurate  wind  data  are  particularly  important 
under  these  conditions. 

Comparison  of  calculated  values  of  the  constituents  of  the  energy  budget  at  the  earth's 
surface  with  measured  values  for  both  the  Davis  and  O'Neill  calculations  showed  the  terres- 
trial and  solar  radiation  calculations  to  be  very  accurate.   Satisfactory  results  were  ob- 
tained for  the  water  vapor  calculation  by  either  specifying  the  latent  heat  flux  or  the 
specific  humidity  at  the  earth's  surface  as  functions  of  time. 

Advection  is  shown  to  be  an  important  energy  transfer  process  and  should  not  be  ne- 
glected in  atmospheric  calculations.   Advection  has  a  significant  effect  upon  the  convective 
heat  transfer  near  the  earth's  surface  even  though  the  rate  of  change  of  temperature  due  to 
advection  may  be  small  in  this  region. 
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1.0   INTRODUCTION 
1 . 1  Motivation 

An  understanding  of  the  energy  transfer  processes  in  the  lower  atmosphere  is  essential 
if   weather    forecasting  is  to  be  accomplished  by  numerical  simulation.   In  order  for  a 
numerical  model  to  be  accurate  in  its  predictions  of  the  weather,  it  must  account  for  the 
various  energy  transfer  processes  accurately.   Accurate  prediction  of  temperature  and 
moisture  fields,  therefore,  requires  accurate  prediction  of  each  of  the  constituents  of  the 
energy  budget. 

Conditions  in  the  lower  atmosphere  are  controlled  by  the  energy  and  water  vapor  bal- 
ances, which  involve  the  complex  energetic  interactions  that  exist  between  the  earth's 
surface  and  the  atmosphere  and  interactions  between  the  large  scale  weather  patterns  of  the 
upper  atmosphere  and  the  smaller  scale  weather  patterns  of  the  lower  atmosphere.   The  pro- 
cesses involved  include  convection,  advection,  solar  radiation,  terrestrial  radiation,  eva- 
poration, condensation,  and  conduction  of  heat  in  the  soil. 

The  predictions  of  both  the  temperature  and  moisture  fields  must  be  accurate  if  the 
numerical  model  is  to  be  capable  of  predicting  the  formation,  dispersion,  and  dissipation 
of  fog  or  stratus.   Because  of  the  relationship  between  vapor  flux  and  latent  heat  transfer, 
the  energy  and  water  vapor  balances  are  intimately  related,  thus  complicating  the  problem. 

There  have  been  relatively  few  attempts  at  quantitatively  predicting  the  temperature 
and  moisture  in  the  lower  atmosphere  as  functions  of  time  and  height.  Generally,  a  diffusion 
approximation  is  used  in  the  numerical  model  to  account  for  the  effects  of  convection,  and 
the  following  numerical  models  use  such  an  approximation. 

Fisher  and  Caplan  (1963)  studied  the  feasibility  of  a  dynamical  model  for  forecasting 
the  development  and  dissipation  of  fog  and  stratus.  The  following  factors  were  considered 
to  be  most  important  in  affecting  the  temperature,  water  vapor  content,  and  liquid  water 
content  of  a  vertical  column  of  air:  (a)  the  vertical  eddy  diffusion  of  these  properties; 
(b)  the  horizontal  advection;  and  (c)  the  latent  heat  associated  with  the  condensation  and 
evaporation.  The  effects  of  radiative  flux  divergence  and  the  evaporation  of  falling  rain 
were  not  included  in  the  model . 

The  temperature,  specific  humidity,  and  liquid  water  content  of  a  vertical  column  of 
the  atmosphere  were  assumed  to  obey  the  following  equation: 


8Q  _  9_ 
9t  ~  3z 


K  f  8Q  1 


3Q    9Q   c  mi 

3x     3y 


where  Q  is  the  atmospheric  variable  being  considered,  K_  is  the  coefficient  of  eddy  dif- 
fusivity,  u  is  the  horizontal  wind  in  the  x  direction,  v  is  the  horizontal  wind  in  the  y 
direction,  and  S  includes  sources  and  sinks  of  the  property  Q. 

An  expanding  vertical  grid  was  employed  in  the  model  in  which  each  grid  length  was  1.25 
times  the  size  of  the  one  below  it,  providing  greatest  resolution  in  the  lowest  layers. 
The  minimum  vertical  grid  length  was  two  meters,  and  the  total  depth  of  the  layer  was  1078 
meters.   The  temperature  at  the  earth's  surface  was  specified,  thus  eliminating  the  need  for 
an  energy  balance  calculation  at  the  earth's  surface. 


The  following  cases  were  considered: 

a)  A  one-dimensional  model  with  the  exchange  coefficient  varying  only  with  height 
(not  time) . 

b)  A  two-dimensional  model  with  the  exchange  coefficients  varying  with  height  only. 

c)  A  one-dimensional  model  with  the  exchange  coefficients  functions  of  the  lapse 
rate  and  height. 

d)  A  two-dimensional  model  with  the  exchange  coefficients  functions  of  the  lapse  rate, 
The  results  were  qualitative  in  nature,  and  the  computed  temperature  and  moisture 

profiles  were  not  compared  with  observational  data.   The  model  demonstrated  the  feasibility 
of  a  dynamical  model  for  forecasting  the  development  and  dissipation  of  fog  and  stratus, 
but  several  improvements  would  have  to  be  made  in  the  existing  model  before  this  could  be 
done  accurately. 

Estoque  (1963)  used  the  following  set  of  equations  to  calculate  the  temperature  and 
moisture  fields: 


3u_  _       1  3_P   _9_  ,„   3u, 

3t     V  ~  p  3x"  +  3x 


■fv-r^^^g  (2) 


(3) 
(4) 


3v     - 
3T=  "  f  u  " 

1_  3P    3     3v 
p  3y    3z  K^i   dzJ 

r      3T    3 
P  p  3t  ~  3z 

3T 

pCp  KH  (  37  +  F  ' 

3R   . 
net 

3z 

3t    3z  L  ^V 

3£ 

3z  J 

(5) 


where  u  is  the  horizontal  wind  in  the  x  direction,  v  is  the  horizontal  wind  in  the  y  dir- 
ection, f  is  the  Coriolis  force,  P  is  pressure,  q  is  specific  humidity,  and  1C.,  K   and  K 
are  the  eddy  diffusivities  for  momentum,  heat,  and  water  vapor  transfers.  R    is  the  net 
longwave  radiation,  and  it  is  taken  as  positive  when  directed  away  from  thesurface.   Only 
the  effect  of  longwave  radiation  is  included  in  Equation  (3)  since  it  is  much  more  signif- 
icant than  the  shortwave  radiation  in  the  lower  atmosphere.   Implicit  in  the  above  equations 
is  the  assumption  of  horizontal  homogeneity  in  the  wind,  temperature,  and  moisture  fields. 
The  pressure  gradients  are  determined  from  synoptic  data  and  are  set  as  constants  in  the 
equations.   The  mean  vertical  motion  is  assumed  equal  to  zero.   Below  50  meters  Estoque 
assumed  that  the  fluxes  of  momentum,  sensible  heat,  and  water  vapor  were  constant  with 
height  and  that  the  turbulent  transfer  processes  were  the  same  for  all  of  the  fluxes.  A 
constant  flux  layer  implies  good  turbulent  mixing  by  the  wind  and  conditions  of  near 
neutral  stability.   A  diffusion  equation  was  used  for  the  soil  temperature  changes.   Steady 
state  was  assumed  at  the  outer  boundaries  of  the  model,  and  the  surface  energy  balance  was 
the  boundary  condition  at  the  soil -air  interface. 

Within  the  lower  50  meters  a  stability  modified  eddy  diffusivity  from  the  log  law  was 
used  for  the  forced  convection  regime,  and  the  eddy  diffusivity  and  profile  shape  result- 
ing from  a  free  convection  analysis  of  turbulent  transfers  was  used  for  the  free  convection 
regime.   Above  50  meters  the  turbulent  transfer  coefficients  were  assumed  to  decrease 


linearly  with  height  to  zero  at  2050  meters.   Values  for  T  and  q  were  computed  between  50 
and  2050  meters  above  the  soil,  then  a  log  law  or(  free  convection  profile  was  fitted  be- 
tween the  surface  and  50  meters  according  to  a  bulk  stability  for  this  layer.   Thus,  there 
was  no  direct  solution  within  the  lower  50  meters. 

Estoque's  computations  resulted  in  fairly  good  agreement  with  surface  energy  balance 
measurements  at  O'Neill  (Lettau  and  Davidson,  1957)  and  with  above  surface  temperature 
fields.   However,  the  agreement  with  the  moisture  fields  was  not  as  good.   Estoque's  assump- 
tion of  a  constant  flux  layer  was  not  consistent  with  the  fact  that  the  temperature  and 
water  vapor  density  changed  with  time  in  this  layer. 

The  divergence  of  longwave  radiation  was  computed  only  above  50  meters,  resulting  in 
negligible  temperature  changes.   However,  radiative  processes  may  cause  significant  temper- 
ature changes  in  the  lower  few  meters  of  the  atmosphere.   The  O'Neill  data  were  obtained 
under  fairly  strong  winds,  so  the  stability  conditions  were  nearly  neutral.   Consequently, 
the  log  law  would  be  a  reasonable  approximation  to  the  profiles  and  transfer  coefficients. 

Zdunkowski,  Henderson,  and  Hales  (1967)  emphasized  the  prediction  of  nocturnal  tem- 
perature changes  in  the  first  few  meters  of  the  atmosphere.   Particular  attention  was 
given  to  the  effects  of  radiative  flux  divergence  and  atmospheric  turbulence  in  this  region, 
but  no  attempt  was  made  to  deal  specifically  with  the  transition  layer.   The  following  equa- 
tions were  used  assuming  horizontal  homogeneity  and  zero  values  for  the  mean  vertical 
velocity  and  local  pressure  change: 
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where  6  is  the  potential  temperature,  T  is  the  atmospheric  temperature,  T  is  the  soil 
temperature,  q  is  the  specific  humidity,  F   is  the  net  radiative  flux,  K    is  the  eddy 
diffusivity  for  water  vapor,  and  K  and  K  are  eddy  diffusivities  for  heat  in  the  atmosphere 
and  in  the  soil.   The  equations  of  motion  were  not  included  in  the  solution. 

The  computed  temperature  profiles  agreed  well  with  the  O'Neill  data  for  the  first  few 
meters  of  the  atmosphere.   Results  indicated  that  the  radiative  flux  divergence  of  the  air 
should  not  be  neglected,  and  the  soil  parameters  must  be  known  to  a  high  degree  of  accuracy. 
Moisture  field  data  were  not  presented  for  comparison  with  the  O'Neill  data. 


Several  improvements  must  be  made  in  existing  numerical  models  before  they  will  be 
capable  of  predicting  temperature  and  moisture  fields  accurately.   However,  once  this  is 
accomplished,  the  model  could  also  be  used  for  studying  the  energy  transfer  processes  in 
the  atmosphere.   The  model  might  also  be  used  for  testing  various  methods  of  climate  control 
which  might  be  used  for  the  removal  of  fog  or  smog. 

Previous  models  have  generally  been  characterized  by  their  very  specialized  nature. 


Estoque's  model,  for  example,  is  intended  for  situations  where  there  is  a  strong  wind,  as 
was  the  case  at  O'Neill,  but  presumably  it  would  not  be  as  accurate  for  situations  with  very 
light  wind.   A  numerical  model  which  is  going  to  be  able  to  forecast  the  weather  accurately 
must  be  capable  of  handling  many  types  of  weather  situations. 

1 . 2  Scope  of  the  Present  Work 

The  primary  goal  of  the  present  work  is  the  development  of  a  numerical  model  which  is 
more  versatile  and  as  accurate  or  more  accurate  than  existing  models  in  predicting  temper- 
ature and  moisture  fields  in  the  lower  atmosphere.   This  is  a  necessary  step  in  achieving  an 
ultimate  goal  of  accurate  weather  forecasting  by  numerical  simulation. 

The  effect  of  horizontal  advection  is  often  ignored  in  numerical  models  on  the  assump- 
tion that  it  is  a  negligible  effect.   A  goal  of  the  present  work  is  the  determination  of 
the  significance  of  advection  on  the  temperature  profile  calculation. 

In  the  present  model,  improvements  are  made  in  the  terrestrial  radiation  calculation, 
and  improvements  are  made  in  the  models  for  the  eddy  diffusion  coefficients  for  heat  and 
water  vapor.   Various  boundary  conditions  aie  used  at  the  earth's  surface  for  the  water 
vapor  calculation  in  an  effort  to  improve  the  accuracy  of  the  latent  heat  transfer  calcula- 
tion at  the  earth's  surface. 

A  diffusion  approximation  is  used  in  the  numerical  model  to  describe  the  effects  of 
convection  and  turbulent  mixing  in  the  atmosphere.   Since  the  wind  field  affects  the  tur- 
bulent mixing,  the  eddy  diffusion  coefficient  depends  to  a  certain  extent  upon  the  wind 
field.   The  wind  field,  therefore,  must  also  be  included  in  the  model.   The  wind  field  may 
either  be  calculated,  as  in  Estoque's  model,  or  it  may  be  specified.   In  the  present  model 
the  wind  field  is  specified.   Consequently,  the  momentum  equations  do  not  have  to  be  solved, 
thus  shortening  the  required  computing  time.   The  wind  field  is  influenced  by  local  terrain 
and  by  conditions  outside  the  region  covered  by  the  model,  so  greater  accuracy  may  be 
achieved  by  specifying  the  wind  field  rather  than  calculating  it.   The  disadvantage  with 
this  method  is  that  wind  data  must  be  available  in  advance  to  cover  the  duration  of  the  cal- 
culation. The  model  may  still  be  used  for  forecasting  purposes  by  using  a  "typical"  set  of 
wind  data  for  the  calculation. 

The  region  of  primary  interest  is  the  first  2,000  meters  above  the  earth's  surface, 
which  is  the  region  where  the  most  rapid  changes  in  temperature  and  moisture  content  occur. 
The  model  covers  a  region  extending  up  to  8,000  meters  above  the  earth's  surface  in  order 
to  include  the  water  vapor  in  this  region  in  the  longwave  radiation  calculation.   The  model 
also  includes  a  layer  of  soil  one  meter  thick,  which  is  needed  for  the  energy  balance  cal- 
culation at  the  earth's  surface. 

2.0   THE  MATHEMATICAL  MODEL 

2- 1   The  Fundamental  Equations 

The  following  equations  are  used  in  the  calculation  assuming  horizontal  homogeneity  and 
negligible  vertical  motion.   The  energy  equation  will  be  modified  later  to  include  an  ad- 
vection term  and  modifications  resulting  from  a  redefinition  of  the  eddy  diffusion  coefficient 


These  modifications  will  be  discussed  in  detail  in  a  later  section. 
The  energy  equation  for  the  atmosphere  is 
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where         9  is  the  potential  temperature 

Ku  is  the  eddy  diffusion  coefficient  of  heat 
H 

R    is  the  net  upward  longwave  radiation 

R   .    is  the  net  downward  solar  radiation 
solar 

C  is  the  specific  heat  of  moist  air  at  constant  pressure 

p  is  the  density  of  moist  air 

z  is  the  vertical  coordinate  measured  positive  upward  from  the  earth's  surface 

t  is  the  time 

P  is  the  pressure 
P  is  a  reference  pressure  (1  atm) 

k  is  the  ratio  (R/C  ) 

R  is  the  gas  constant  for  moist  air. 
The  potential  temperature  is  defined  by 

P 
9  =  T  (  p2.)K  (10) 

where  T  is  the  temperature  of  the  air.  The  potential  temperature  is  the  temperature  which 
would  result  if  the  air  were  brought  isentropically  to  a  pressure  of  one  atmosphere. 
The  diffusion  equation  for  water  vapor  is 


3^   _3_       3^         lsr<*M«rv.&         ?J*.     ^      /*W^ 

3t  ~  3z  [  p  \   3z  J  l   J 


P 


where  q  is  the  specific  humidity  (gm/gm  air)  and  1C.  is  the  eddy  diffusion  coefficient  of 
water  vapor. 

The  energy  equation  for  the  soil  is 

3T  3T 
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g        g 

where  T  is  the  temperature  of  the  soil,  K  is  the  thermal  conductivity  of  the  soil,  and  z 

g  F  g  g 

is  the  vertical  coordinate  for  the  soil  measured  positive  downward  from  the  earth's  surface. 

The  effect  of  precipitation  and  evaporation  of  liquid  water  is  not  included  in  these 

equations.   Each  time  new  values  are  calculated  for  temperature  and  specific  humidity,  a 

test  is  made  to  see  if  the  air  is  supersaturated,  and  the  temperature  and  specific  humidity 

are  adjusted  accordingly. 

It  is  assumed  that  the  heat  capacity  of  the  soil  C   ,  the  density  of  the  soil  p  ,  and 

vi  gr»  g» 

the  thermal  conductivity  of  the  soil  K  are  constant.   It  is  also  assumed  that  air  behaves 

g 
as  an  ideal  gas,  thus  the  equation  of  state  for  air  is 
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P  =  pRT  (13) 

The  gas  constant  for  moist  air  is  variable  because  the  specific  humidity  of  water  vapor 
in  the  atmosphere  is  variable.  The  gas  constant  R  can  be  expressed  in  terms  of  the  gas 
constant  for  dry  air  R,  and  the  gas  constant  for  water  vapor  R   by  the  relationship 
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where  q  is  the  specific  humidity  of  water  vapor.   Assuming  that  R,  and  R   are  constant, 
we  get 


(1  +  0.61q)  Rd 


(15) 


The  heat  capacity  of  moist  air  at  constant  pressure  is  also  variable  because  of  var- 
iations in  the  specific  humidity  of  water  vapor.   The  heat  capacity  of  moist  air  C  can  be 
expressed  in  terms  of  the  heat  capacity  for  dry  air  C  ,  and  the  heat  capacity  of  water 
vapor  C    by  a  relationship  similar  to  that  for  the  gas  constant: 


l  +  (; 


pwv 


i)  q 


>d 


(16) 


Assuming  that  C  ,  and  C^irir  are  constant,  we  get 


pa      pwv 

C   =  (1  +  0.90q)  C 


Pd 


(17) 


2. 2  Boundary  Conditions 

At  the  upper  boundary  of  the  atmospheric  region,  the  heat  flux  and  water  vapor  flux  are 
set  equal  to  zero.   Expressing  this  in  terms  of  the  variables,  we  gat 
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At  the  lower  boundary  of  the  soil  layer,  the  heat  flux  is  set  equal  to  zero,  which  means 

3T 
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At  the  earth's  surface  the  temperature  is  calculated  using  an  energy  balance,  and  the 
specific  humidity  is  calculated  using  one  of  several  methods.   Among  the  various  boundary 
conditions  used  in  the  water  vapor  calculation  are  the  following: 

1)  the  specific  humidity  is  held  constant 

2)  the  relative  humidity  is  held  constant 

3)  the  specific  humidity  is  specified  but  allowed  to  vary  with  time 

4)  the  latent  heat  flux  is  specified  and  allowed  to  vary  with  time 
These  methods  are  used  as  a  preliminary  step  before  trying  more  sophisticated  methods 


of  calculating  the  specific  humidity  at  the  earth's  surface.   The  specific  humidity  may  be 
calculated  using  a  moisture  balance  at  the  earth's  surface,  and  more  advanced  models  may  in- 
clude the  effects  of  plant  physiological  processes.   Processes  such  as  stomatal  closure  may 
significantly  affect  the  transpiration  of  water  vapor  from  plants,  thereby  influencing  the 
specific  humidity  at  the  earth's  surface. 


2. 3  The  Longwave  Radiation  Calculation 

Define  the  optical  thickness  for  water  vapor. 
z„ 


pq  dz 


(20) 


'1 


where  z_  is  the  height  at  the  top  of  the  atmospheric  region,  and  z   is  the  height  correspond- 
ing to  this  value  of  w.   The  optical  thickness  represents  the  amount  of  precipitable  water 
contained  between  the  top  of  the  atmospheric  region  and  the  height  z..  .   So  to  is  zero  at  the 
top  of  the  atmospheric  region  and  increases  moving  downward. 
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Figure  1.   Longwave  radiation  from  an  elemental  layer  of  the 
atmosphere. 

The  radiative  flux  dF  of  an  infinitely  extending  horizontal  elemental  layer  of  optical 
thickness  dco  and  temperature  T  arriving  at  a  horizontal  receiver  after  traversing  the  ab- 
sorbing mass  (V  between  the  detector  and  the  emitting  layer  is 
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where  \>   is  the  frequency,  k   is  the  absorption  coefficient,  and  B  =  1.66  =  sec 
B   is  the  black  body  radiation  for  a  frequency  interval  dv . 

A 


B 


B  dv  =  aT 

v 


(22) 


where  a  is  the  Stefan-Boltzmann  constant 


Since  the  temperature  varies  with  height,  and  there  is  a  value  of  to  associated  with 
each  value  of  z,  the  temperature  is  in  effect  a  function  of  a>.   Consequently,  B  is  also  a 
function  of  to. 

Make  the  following  definition  for  the  transmissivity  t. 
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and 
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We  are  interested  in  the  net  upward  radiative  flux  at  any  reference  height.  Once  that 
is  known,  the  rate  of  change  of  temperature  due  to  the  absorption  of  radiation  can  be  cal- 
culated by  taking  the  divergence  of  the  net  radiation  at  that  height.   Let  w  be  the  optical 
thickness  at  the  reference  height  and  let  W  be  the  optical  thickness  at  the  earth's  surface. 
(Refer  to  Figure  2.) 


reference 
height 


"  to 


earth's  surface 


Figure  2.  The  optical  thickness  at  various  heights  above  the 
earth's  surface. 

Using  the  definition  of  t  given  above,  the  radiation  which  is  emitted  by  the  earth's 
surface  and  is  transmitted  through  the  intervening  atmosphere  to  the  reference  height  is 


e   B  t(W  -co  ) 
s  s      r 


where  c   is  the  emissivity  of  the  earth's  surface  and  B   is  the  black  body  radiation  cor- 
responding to  the  temperature  of  the  earth's  surface. 

The  radiation  emitted  by  the  air  layer  below  the  reference  height  is  found  by  inte- 
grating equation  (21)  over  that  layer. 
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Similarly,  the  radiation  emitted  by  the  air  layer  above  the  reference  height  is 
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The  radiation  which  is  emitted  by  the  entire  atmospheric  region,  reflected  by  the  earth's 
surface,  and  reaches  the  reference  height  is 
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The  net  upward  radiative  flux  at  the  reference  height  is  the  sum  of  these  fluxes, 
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Integrating  by  parts  and  substituting  dB(w)  =  -^  dw,  we  get 
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The  divergence  of  the  net  radiation  is 
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At    the  reference  height, 
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Noting  that  dB(co)  =  —  dco,    we  get  the  following  expressions 
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In  order  to  carry  this  development  further,  it  is  necessary  to  specify  a  form  for  x  as 
a  function  of  w.   Crawford  (1965b)fit  a  curve  to  the  data  of  F.  A.  Brooks  (1941).   These 
data  include  the  effect  of  CO.  concentrations  of  about  0.03  percent.   That  same  expression 
for  x  is  used  here. 
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Table  1:   Coefficients  for  transmissivity  of  terrestrial  radiation 

k  a, 


k 

0.01900 
0.07000 
0.18000 
0.26400 
0.17445 
0.29255 


Mk 

6000.00 

1700.00 

270.00 

30.00 

2.50 

0.15 
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The  factor  of  1.66  for  slab  emission  is  included  in  6,  .   Substituting  this  expression 
into  the  equation,  we  get 
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2.4  The  Solar  Radiation  Calculation 
Let  I  =  solar  constant 
<t>  =  latitude 
6  =  solar  declination 
h  =  hour  angle 
u  =  turbidity  factor 
The  turbidity  factor  accounts  for  absorption  and  reflection  of  solar  radiation  in  the 
upper  atmosphere.   Then  the  solar  radiation  incident  at  a  horizontal  receiver  at  the  top  of 
the  atmospheric  region  is 


R   ,    =  yl„  (sind)  sin6  +  cost))  cos6  cos  h)  . 
solar     0  v 

If  the  zenith  angle  of  the  sun  is  0,  then  R  ,    is  also  given  by 
&  solar        b 

R   ,    =  yl^  cos  0. 
solar     0 


(32) 


(33) 


Comparing  equations  (32)  and  (33) ,  we  see  that 

cos  0  =  sin  <j>  sin  6  +  cos  cf>  cos  6  cos  h.  (34) 

If  t   ,    Too  sec  0)  is  the  transmissivity  of  water  vapor  for  solar  radiation,  then  the 
solar  j  r 

solar  radiation  which  reaches  a  horizontal  receiver  after  traversing  an  absorbing  mass  w  is 


R   -,    =  ul„  cos  0  t   n   (no  sec  0)  . 
solar     0        solar 


(35) 


An  analytic  expression  can  be  found  for  t   ,   by  fitting  a  curve  to  the  data  by  H.H.  Kimball 
1  r  solar  '       & 


(1931).   An  expression  of  the  following  form  was  used  for  t   , 


ar 


t  „lQJu  sec  0)  =  a^  +  I       a  exp(-b,  u>  sec  0) 


solar 


k-1  k 


(36) 
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Table  2:   Coefficients  for  transmissivity  of  solar  radiation 

k               ak  Sk 

1  0.0585780  0.10 

2  0.0300355  0.50 

3  0.0164936  1.00 

4  0.0436047  5.00 

5  0.0094854  10.00 
a„  0.8418028 


*0 

The  divergence  of  the  solar  radiation  is  given  by 

solar  _    solar  3co 
"3z        3co      3z 


From  equation  (35),  we  see  that 


(37) 


,  SOlar  =  pi  cos  0  /solar   (co  sec  0)  (38a) 

3co  o       3co 


or 


solar     T   ^ 


-yl   E   a,  b,  exp  (-b,  co  sec  0)  (38b) 


3co      ~~   M  o  ,  ,  "k  "k   ^  v  k 
k=l 

2 . 5  The  Eddy  Diffusion  Coefficients  • 

T.  V.  Crawford  (1965)  developed  an  expression  for  the  eddy  diffusion  coefficient  which 
is  fitted  to  heat  flux  data  and  assumes  a  stability  modified  log-law  velocity  profile.   His 
expression  for  K„  has  the  advantage  that  it  allows  a  smooth  transition  between  the  forced 
convection  regime  and  the  free  convection  regime.   It  is  assumed  that  the  eddy  diffusion 
coefficient  for  water  vapor  equals  the  eddy  diffusion  coefficient  for  heat,  which  is  a 
reasonable  assumption  based  upon  measurements  of  K  and  K„  near  the  earth's  surface. 

The  Richardson  number  is  used  as  a  measure  of  thermal  stability,  and  it  is  defined  here 
as 

(■5—  +  n 

R,  =  I- (39) 

8U  2 
^zJ 

where  g  is  the  acceleration  of  gravity,  r  is  the  adiabatic  lapse  rate  (0.01°C  m  ),  and  U 
is  the  horizontal  wind  velocity. 
The  heat  flux  is  given  by 


and 


"  =  "p  Cp  KH  <!r  +  r  ^  (40) 


2      2 
k   (z-d) 


KH  = 


3U_ 

3z 


(41) 


<j>(R.) 
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where  k  is  von  Karman ' s  constant  and  d  is  the  displacement  height.   The  function  4>(R.)  is  a 
stability  correction  to  the  eddy  diffusion  coefficient  for  neutral  conditions.   Crawford  ob- 
tained expressions  for  4>(R.)  for  the  forced  and  free  convection  regimes  by  fitting  a  curve 
through  heat  flux  data  for  the  first  10  meters  of  the  atmosphere.   For  the  free  convection 
regime  (R.<0),  k  =  0.426  and 


>(R.)  =  (1-0.25  R.) 


-0.60 


(42) 


For  the  forced  convection  regime  (R.>0),  k  =  0.417  and 


>(R.)  =  (1+0.24  R.) 
l    v        \J 


0.40 


(43) 


This  expression  for  K  is  only  good  near  the  earth's  surface  where  the  log-law  velocity 
profile  is  a  good  approximation,  so  another  expression  is  needed  for  the  region  away  from 
the  earth's  surface. 

J.  W.  Deardorff  (1967)  suggests  the  following  definition  of  K  for  the  region  above 
50  meters . 


H  = 


-  Cp  KH  <H 


V 


(44) 


where  y   is  a  small  positive  quantity  which  allows  a  counter-gradient  upward  heat  flow.   It 
has  been  noted  by  Deardorff  that  a  counter-gradient  upward  flux  appears  to  occur  above  the 


lowest  100  m  or  so  when  0  <  - —  <  0.65  x  10~~ 

3z 

The  formulation  adopted  by  Deardorff  is 


'C  cm   in  clear  air. 


KH  =  KH  /  C1  +  Cl  RP  £°r  H  "  \ 


(45) 


K^  =  KH  +  C2{  1-exp 


a(97  *  V 


i  ,   36  < 
J  £or  3T  "  \ 


(46) 


R'.  is  an  adjusted  Richardson  number  given  by 


R'.  = 

i 


,36     » 
;(37  '  Yc} 


,3u.2    ^3v~.2 

^    +  (3r} 


(47a) 


where  u  is  the  horizontal  wind  in  the  north-south  direction  and  v  is  the  horizontal  wind  in 

* 
the  east-west  direction.   K   C.  and  C,  are  constants  to  be  determined  by  comparison  of  cal- 

H    1       2 

culated  temperature  profiles  with  observed  ones.   For  zero  turbulent  heat  flux,  it  is  seen 
that  ■£  =  y  and  K'  =  K*   In  order  that  3K'  /  3(36/3z)  be  continuous  at  r-=  y  ,  the  constant 

oZ     C        H     H  H  a  Z     C 

a  must  be  given  by 
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^1 


<lr)2 


(47b) 


a  9 
where  the  subscript  N  indicates  evaluation  at  the  "no-flux"  level  where  -r—  =  y    , 

such  levels  occasionally  occurred,  the  upper  one  was  selected  to  evaluate  a. 

* 
Deardorff  arrived  at  the  following  values  for  K   C.,  and  C„: 

*         5   2 
K„  =  6  x  10  cm  /sec 

H 


'1 


50 


When  two 


f\  9 

C?  =  6  x  10  cm  /sec. 
It  is  possible  to  define  an  eddy  diffusion  coefficient  which  will  be  good  for  all 
heights  by  combining  the  two  formulations  above.   The  definitions  of  K'  and  R.'  can  be  ex- 
tended to  include  the  lowest  50  meters  of  the  atmosphere  by  setting  y     equal  to  zero  in  this 
region.   When  this  is  done,  K'  is  defined  in  the  first  50  meters  of  the  atmosphere  by 


»=-'VJ5 


(48) 


The  gradient  of  potential  temperature  is  related  to  the  gradient  of  absolute  temperature  by 
the  expression 

(49) 


96_ 
3z 


■£)' 


£♦' 


By  substituting  this  expression  into  equation  (48)  and  comparing  it  with  equation  (40),  we 
see  that 


K'h  =  (P~)   KH 
o 


(50) 


A  relationship  between  R.  and  R!  is  found  by  substituting  equation  (49)  into  equation 
(47)  with  y     equaling  zero  and  comparing  the  result  with  equation  (39).   When  this  is  done, 
we  see  that 


Ri =  <H 

o 


,3u..2    ,3v.2 


(51a) 


3U  2 
l9zJ 


p   K  p   K 

For  the  region  below  50  meters  (p— )   =  1.   By  assuming  (■=— )   =  1 ,  we  get  an  expression 

for  the  eddy  diffusion  coefficient  near  the  earth's  surface.   Let  K'  represent  the  eddy 

1 
diffusion  coefficient  for  this  region,  then 


where 


k2(z-d)2 
\=        *(R  ) 


8U_ 

3z 


(51b) 


,,3u.2   ,-9v.2 
R.  =  Jl i^_  R; 

l3zJ 


(51c) 
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Let  K'   represent  the  eddy  diffusion  coefficient  for  the  region  above  50  meters,  then 
KH2  =  KH  /  <1+C1  RP  for  H  >\  (51d) 


and 

K^  =  K*  +  C2{1  -  exp  - 


,36     . 
a(3T  "  V 


l  ^   39 

}  f or  j-  <  yr  (5le) 


The  eddy  diffusion  coefficient  is  defined  by  K'   near  the  earth's  surface  and  by  K' 

1  2 

above  a  height  of  about  50  meters.   The  eddy  diffusion  coefficient  can  make  a  smooth 

transition  between  K'   and  K'   by  using  a  combination  of  the  two  definitions  in  a  transit- 

1        2 
ion  region.   If  the  transition  occurs  between  the  heights  z  and  z  ,  where  z  <  z    ,    then  K' 

is  defined  by 

K'   =  K '     for  0  <  z  <  z,  with  y     =   0 
H    H  -    -   1       'c 


with  >   =0  where 


KH  "  h    \   +  h   \  f0r  Zl  <   Z  <  Z2 


Z--Z  z-z 

C,  =  and     Ct  =  

1   z2-zx  ^2   z2-Zl 


K'h  =  K'h2     f0r  Z  >  Z2 

The  heights  z  and  z  are  chosen  to  be  20  meters  and  50  meters,  respectively.   The 
height  of  50  meters  is  chosen  because  the  earth's  surface  has  little  effect  on  the  turbulent 
mixing  at  this  height,  and  K'   applies  to  regions  outside  the  influence  of  a  physical  boundary, 

A  log- law  velocity  profile  was  assumed  in  the  development  of  K'  ,  which  applies  to  the 

1 
region  near  the  earth's  surface.  The  height  of  20  meters  is  chosen  as  an  estimate  of  the 

maximum  height  at  which  the  earth's  surface  significantly  affects  the  turbulent  mixing  and 
the  wind  profile  can  still  be  approximated  by  a  log-law  velocity  profile.   The  heights  z 
and  z  are  not  strictly  defined,  and  they  may  be  varied  during  the  calculation  in  future 
models  according  to  the  thermal  stability  and  wind  conditions. 

There  are  several  ways  in  which  the  two  definitions  of  K'   can  be  fitted  together.   In- 
stead of  combining  the  two  definitions  over  a  transition  region  between  z  and  z  ,  as  was  done 
above,  the  two  definitions  could  be  fitted  together  at  a  certain  given  height.   If  this  were 

done.  K'   and  K'   would  be  forced  to  have  the  same  values  at  that  point  along  with  the 

1        2 
values  of  their  first  derivatives  with  respect  to  height.   By  doing  this,  the  values  of  K' 

1 

and  K'    would  have  to  be  changed  in  the  vicinity  of  the  point  where  they  were  fitted. 

Rather  than  do  this,  the  method  used  above  allows  K'   and  K'u  to  remain  unchanged,  and  it 

i        o 
enables  a  gradual  transition  from  an  eddy  diffusion  coefficient  which  is  influenced  by  a 

boundary  to  one  which  is  not  influenced  by  a  boundary.   In  the  atmosphere,  the  effect  of  the 

earth's  surface  on  turbulent  mixing  becomes  less  as  you  get  farther  away  from  the  surface. 
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This  transition  is  gradual  and  not  abrupt,  which  is  another  reason  for  using  the  method  of 
weighting  factors. 

In  the  original  energy  equation,  the  heat  flux  was  represented  by 

H  =  -d  C   KH  |i 
p   H  9z 

The  energy  equation  must  now  be  modified  because  K'  is  defined  by 
H  =  "p  Cp  KH  C||  -  Yc) 


3  fi  3  fi 

This  means  that  K'  (^ y  )  must  be  substituted  for  K „(•»— )  in  the  energy  equation.   When 

H   a  Z      C  H  o  Z 


this  is  done,  we  have 


0  cp  at   9z 


0  Cp  KH  <i  -  V 


P   K   3R   . 
,  o.   ,  solar 

+  *r>    {Tz 


3R 


net 


3z 


(52) 


2.6  Advection 

The  energy  equation  as  developed  up  to  this  point  in  the  presentation  contains  no  ad- 
vection terms.   Assuming  negligible  vertical  motion,  this  study  is  restricted  to  horizontal 
advection.   The  presence  of  an  horizontal  temperature  gradient  affects  heat  transport  by 
turbulent  mixing  as  well  as  by  advection.   Since  the  horizontal  temperature  gradients  studied 
in  this  work  are  less  than  0.5°C  km   ,  its  effect  on  turbulent  mixing  is  negligible.   When 
advection  is  included,  the  energy  equation  becomes: 


„   ,-,.36    30,.    3     „   ,.,,38 
p    3x   3t     3z     p  H  v3z 


^ 


f       0-! 


36 


3R 


solar 


3R 


net 


3z 


3z 


(53) 


where  U  is  the  magnitude  of  the  horizontal  wind  vector  and  —  is  the  gradient  of  potential 

oX 

temperature  in  the  direction  of  the  wind  vector. 


3.0  COMPUTATIONAL  METHODS 

3 . 1   The  Grid  in  the  Atmospheric  Region 

An  expanding  vertical  grid  is  used  in  the  atmospheric  region  in  which  each  grid  length 
is  1.25  times  the  size  of  the  one  below  it,  thus  providing  greater  resolution  in  the  lower 
layers.   The  atmospheric  region  extends  up  to  a  height  of  8070  meters  and  contains  32  grid 
points.   The  minimum  grid  length  is  2.0  meters,  and  25  grid  points  lie  within  the  lowest  2000 
meters,  which  is  the  region  of  primary  interest.   Let  z(j)  represent  the  height  of  the  j 
grid  point.   Then 

z(l)  =  0.0  meters 

z(2)  =  2.0  meters 

z(3)  =  4.5  meters 

etc . 
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3. 2  The  Grid  in  the  Soil  Region 

An  expanding  vertical  grid  is  also  used  in  the  soil  region  in  which  each  grid  length 
is  1.25  times  the  size  of  the  one  above  it.   The  soil  region  extends  to  a  depth  of  1.08 
meters  and  contains  13  grid  points.   The  minimum  grid  length  is  2.0  centimeters.   A  maximum 
depth  of  1.08  meters  is  used  because  the  diurnal  temperature  changes  are  small  at  this  depth 
compared  to  the  temperature  changes  at  the  surface.   Let  z  (j)  represent  the  depth  of  the  j 
grid  point.   Then 

z  (1)  =  0.0  centimeters 

g 
z  (2)  =  2.0  centimeters 

g 
z  (3)  =  4.5  centimeters 

g 
etc. 

3. 3  Expressing  Derivatives  in  Finite-Difference  Form 

An  expanding  vertical  grid  has  the  advantage  of  improving  resolution  near  the  earth's 
surface,  but  it  presents  a  problem  when  it  comes  to  writing  the  equations  in  finite  dif- 
ference form.   The  usual  method  of  expressing  derivatives  as  the  difference  between  values 
of  a  variable  at  certain  grid  points  divided  by  the  distance  between  those  grid  points  re- 
sults in  expressions  which  are  not  centered  on  grid  points  when  an  expanding  grid  is  used. 
Most  finite-difference  methods  are  not  applicable  to  expanding  grids  without  sacrificing 
accuracy.   It  is  possible  to  get  around  this  problem  by  a  transformation  into  another  co- 
ordinate system  where  the  grid  points  are  evenly  spaced.  To  determine  what  the  trans- 
formation must  be,  we  begin  by  developing  a  functional  form  for  z(j).  The  height  of  the 
j   grid  point  in  the  atmospheric  region  is  given  by 


z(j)  =  800.0 


(1.25)3"1 


-  1 


cm  (54) 


If  we  now  consider  j  to  be  a  continuous  variable  instead  of  a  discrete  variable,  j  is  no 
longer  a  subscript  but  an  independent  variable.   For  every  value  of  z  there  is  a  corres- 
ponding value  of  j .   Grid  points  correspond  to  integer  values  of  j ,  so  the  expanding  grid 
in  "z-space"  has  been  transformed  into  an  evenly  spaced  grid  in  "j -space". 

We  now  develop  an  expression  for  the  partial  derivative  in  finite  difference  form. 
Since  j  is  a  continuous  variable, 


8<j)  _  d$_     dj_ 
9z  ~  3j   3z 


Solving  for  j  using  equation  (54),  we  get 

ln(l  +  z/800.0) 
3  ln(1.25) 


(55a) 


+  1  (55b) 


Taking  the  partial  derivative  with  respect  to  z  and  substituting  equation  (54)  into  the 
result,  we  get 
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ii 

3z 


800    (In   1.25)    (1.25) 


j-l 


(56) 


We  approximate  the  partial  derivative  of  <J>  by 


3<j>   _   A<$> 
3j    "   Aj 


(57) 


Substituting  equations  (56)  and  (57)  into  equation  (55a), we  get 


3^         <l>(j+l)  -  4>(j) 
l3zJj+l/2 


800  (In  1.25)  (1.25) 


j-1/2 


(58) 


We  develop  an  expression  for  the  second  derivative  as  follows 


9  (j)  _  dj_   3 ,3 


(tr) 


3z 


2   3z  3j  '-Sz 


(59a) 


+■  v» 
The  second  derivative  evaluated  at  the  j   grid  point  is  approximated  by 

1 


A), 


r3<K        f3fK 

l3zJj+l/2  "  l-3z''j'l/2 


ann  rin  i  7^\    m  9mJ~ 


800  (In  1.25)  (1.25)J     (j  +  1/2)  -  (j  -  1/2) 


(59b) 


Making  substitutions  for  the  first  derivatives,  we  get 


ti±)      =   <t>(3+l)    -   2.25   <j>    (j)    +   1.25   <|>    (j-l) 
3z2   j  64  x   104    (In  1.25)2    (1 .25)2;j  ~3' 2 


(59c) 


V 


3.4   Solution  of  the  Water  Vapor  Equation 

The  differential  equation  to  be  solved  is: 


pfMr^if) 


(60) 


Let  q  ( j )  =  the  value  of  q  at  the  j   grid  point  and  the  n   time  step 
and    At  =  the  time  step 


i  +  1 


Height 


-1 


._Z 


n-1    n    n  +  1 
Time 

Figure  3:   The  time-height  grid 
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grid  points  used  in 
the  calculation 


Assume  that  q  (j)  is  known  for  all  values  of  j .   We  express  the  differential  equation 
in  finite  difference  form,  and  solve  for  q   (j)  in  terms  of  q  (j-1),  q  (j),  q  (j+1), 
q   (j-1),  and  q    (j+1).   The  points  used  are  shown  in  Figure  3. 

Equation  (60)  is  approximated  by  the  finite-difference  system 


n+1  n,.. 

q      (j)-q  (j)  =  r 

At  ^ 


fIt<"<v& 


n+1 


+    (1-?) 


JD 


19     t    v    M\ 

p"  37   (p    Sf  J? 


(61a) 


J  3 


where   ^   is  a  weighting   factor   in  the   interval   0   <   £   <   1.      If  5=0,    the   scheme   is  purely 
explicit,    and   if  £  =   1,    the   scheme   is  purely   implicit. 
From  equation    (58) 


iFfe&^fc 


n+1 


1 


j  Pn+1(j)    64   x   104    (In   1.25)2    (l^S)1""1 

•{pn+1  cj+i/2)  Kin+1  o*i/2)  qn+1(j+1):C(j) 

(1.25)J      /Z 


n , ..      n 


-pn+1(j-l/2)    lL.n+1    (j-1/2)    S    <JH    \>-)l    }  (61b) 


Assume   1C        (j)    =   K.     (j)    and   p        (j)    =   p    (j).      This  may  be  done   for  small   time   steps 


because   1C    and   p  vary   slowly  with  time.      Then 
13      ,      ,,      3c 


I    d       t      v      dq. 

p    97   (P    *V    37} 


j  pn(j)      64   x   104    (In  1.25)2    CL  25)  ^    * 

•{pn(j+l/2)    Kvn(j+l/2)    qn+1CJ^l)-q"^(J) 


n+1,.-.      n+1,.    . . 

-Pn(j-i/2)  ^(j-i/2)  q    w-q  .  [)f\ 

(1.25)J      ' 


(61c) 


Substituting  expressions  for  the  derivatives  on  the  right  hand  side  of  equation(61a)  and 
multiplying  through  by  At,  we  get 


n+1,..   n 

q   (j)-q  0)  = 


2At 


pn(j)  128  x  104  (In  1.25)2  (1 .25)2;i  "3/2 


{pn(j+l/2)  ^0+1/2) 


_ ,  n+1 , .  , .  n+1 , ... 
C(q   (j+l)-q   0)) 


Let 


+  (1-5)  (qn(j+D-qn(J)) 


-  1.25  pn(j-l/2)  ^"(j-1/2) 


«(j) 


5(qn+1(J)-qn(J-D)  +  d-C)(qn(j)-qn(j-D) 


At 

128  x  104  (In  1.25)2  (1 . 25) 2J "3/2 


(61d) 
(61e) 
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We  arrange  terms  in  the  following  order 

V 


n+1,. 


_2^iii_pn(j  +  1/2)  Kn  (j+1/2)  qn+1  (j+l)    S&l-  2 .50Cpn(j  -1/2)  kJ  (j-1/2)  q-*(j-l) 


n  . 

o    (j) 


P  (j) 


*lJ 


U2i^L-{    pn(j+l/2)  K"  (j+1/2)  +  1.25pn(j-l/2)  k"  (j-1/2) 


n, . 
P  (j) 


n+1 

q   0) 


■-2^p-   (l-4)on  (j+1/2)  Kj  (j+1/2)  q"  (j+l)  +  2.50  ^lil_(l-Opn(j  -1/2)  kJ  (j-1/2) 
P  (j)  P  (j) 


q  U-D 

1-2  Sill.  (l-O  |pn(j  +  l/2)  K"  (j+1/2)  +  1.25pn(j-l/2)  Kj  (j-1/2)} 
P  (j) 


qn(j)  (62) 


This  equation  is  of  the  form 


where 


and 


■\(j)  qn+1  Cj+l)  +  By(j)  qn+1  (j)  .  Cv(j)  qn+1  (j-1)  =  Dy(j) 


Av(j)  =  2^(j)pn(j+l/2)Kj(j+l/2)/pn(j) 


Bv(j)  =  l+25o(j)  j  pn(j+l/2)K^(j+l/2)  +  1.25  pn(j -1/2)k"(j -1/2) 

Cv(j)  =  2.50Ca(j)pn(j-l/2)  Kj  (j -l/2)/Pn(j ) 

r~ 
Dy(j)  =  2(l-Oct(i)!pn(j+l/2)  Kj  (j+l/2)qn(j+l) 

+1.25pn(j-l/2)  K^  (j-l/2)qn(j-l) 


/P"(j) 


/Pn(j) 


+  {l-2(l-Oa(j) 


(63) 


n,». 


pn(j+l/2)  Kj  (j+1/2)  +  1.25pn(j-l/2)  Kj  (j-l/2)|/pn(j)}qn(j) 


■r      1 


With  this  scheme,  stability  is  insured  if  ■=■<£<  1 .   For  this  range  of  £,  a  value  of 

1 

y  provides  the  greatest  accuracy,  so  that  value  was  used  for  £. 

Use  the  algorithm 


qn+1  (j)  =  E.,(j)  q"*1  (j-1)  +  Fu(j) 


(64) 


The  coefficients  Ev(j)  and  F  (j)  are  determined  by  substituting  E  (j  +  l)  q(j)  +  F  (j  +  l)  for 
q(j+l)  and  rewriting  the  equation  in  the  form 


n+1  . 
q   (J) 


sw 


Bv(j)-Av(j)Ev(j+l) 
Comparing  equations  (64)  and  (65),  we  see  that 


cv(j) 


n+ir.  n    vj)*yj)yj+i) 

<*    "-U  +  Bv(j)-Av(j)Ev(j+l) 


LV(J)  =  Bv(j)-Av(j)Ev(j+l) 


(65) 


(66a) 


and 
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Dv(J)^(J)Fv(j+l) 
VUJ  "  Bv(j)-Av(j)Ev(j+l) 


(66b) 


Apply  the  boundary  conditions  on  q  to  equation  (64).   At  the  upper  boundary  -r-*-  =  0. 
This  means  q    (j)  =  q    (j-1),  where  j  is  the  grid  point  at  the  boundary.   From  equation 
(64) ,  we  get 


Ey(j)  =  1.0  and  Fy(j) 


0.0 


(67) 


Note  that  E  (j)  and  Fv(j)  depend  on  the  values  of  the  variables  at  the  present  time 
step  only,  thus  enabling  E  (j)  and  F..(j)  to  be  calculated  without  knowing  q   ( j ) .   The 
value  of  the  specific  humidity  at  the  earth's  surface  must  be  known  for  the  next  time  step 
before  q   (j)  may  be  calculated.   According  to  this  scheme,  E  (j)  and  F..(j)  are  calculated 
moving  down  the  grid  from  top  to  bottom,  then  q   (j)  is  calculated  using  equation  (64) 
moving  up  the  grid. 


3.5  Solution  of  the  Energy  Equation  for  the  Atmosphere 
The  energy  equation  for  the  atmosphere  is 


PCP  (U3x-  +  9T}  =  37 


pCp  KH  <H  -V 


3R   .     3R 
solar    net 


3z 


3z 


P  K 

(— ) 


Dividing  through  by  pC  and  rearranging  terms, 


(68a) 


If  =  w~  &  (pCP  KA  ff}  -  plr  k  (pCP  kh  V 
p  F  p  F 


.         3R      .  3R     ^ 

1  solar  net 

pC        3z 
P  L 


3z 


P       K 

(-2.)     -u™ 
KP  J  3x 


(68b) 


Assume  that  C  is  constant  over  the  range  of  the  six  grid  points  used  in  the  scheme, 


wh 


1   3 


ich  covers  one  time  step.   When  the  term  -  —=—  -r—   (pC  K'y  )  is  expressed  in  finite  difference 

F  pC   3z  tH  p  H'c'      r 

form  using  equation  (58),  we  get  " 


p(j+l/2)  K'H(j+l/2)  Yc(j+l/2)  -  p(j-l/2)  K^  (j-1/2)  Yc(j-l/2) 


p(j)  800  (In  1.25)  (1.25) 


j-l 


Using  a  development  similar  to  that  used  in  the  water  vapor  calculation,  letting  £=  ■=-, 
assuming  o        (j)  =  p  (j)  and  K   (j)  =  K„(j),  and  arranging  terms  in  the  form 


-AH(j)  en+1(j+i)  +  bh(j)  en+1(j)  -  cH(j)  en+1(j-i)  =  dh(j), 


(69) 


we  get  the  following  expressions: 

AH(j)  =  a(j)p(j+l/2)  K^  (j+l/2)/p(j) 


BH(j)  =  1  +  a(j) 


p(j+l/2)  K^  (j+1/2)  +  1.25  p(j-l/2)  K^  (j-1/2) 


/P(j) 
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CH(j)    =    1.25   a(j)p(j-l/2)    K^    (j -l/2)/p (j) 


DH(j)    =   a(j) 


p(j+l/2)    K^    (j+l/2)9Cj+D    +    1-25    pCj-1/2)    K^    (j -1/2)  6  (j -1) 


/P(j) 


^(j){l-a(j) 
i~9R 


>At{ 


9z 


p(j+l/2)    K^    (j+1/2)    +    1.25    p(j-l/2)    K^    (j-1/2) 
p(j)C 


/P(j)} 


9R    . 

solar  net 


9z 


P         K 


HOI    C&). 


p(j+l/2)    K^    CJ+l/2)Yc(j+l/2)    -   P(j-l/2)    K^    (j -1/2)  yc  (j -1/2) 


p(j)    800    (In    1.25)    (1.25) 


j-l 


We  use  the  algorithm 
„n+l 


,n+l 


(j)  =  EH(j)  6"  *(j-l)  +  FH(j) 


(70) 


which  leads  to  the  following  expressions 
EH(j)  = 


cH(j) 


and 


BH(j)-AH(j)EH(j+l) 


FH(j)  =  »HM+WHW 
BH(j)-AH(j)EH(j+l) 


(71a) 


(71b) 


Applying  the  boundary  condition  on  6  at  the  upper  boundary  of  the  atmospheric  region 

9z 

we  get  the  relationship  9    (j)  =  6    (j-l)  where  j  is  the  grid  point  at  the  boundary.   From 

equation  (70)  we  see  that  E„(i)  =1.0  and  Fu(j)  =0.0 

H  n 

As  in  the  water  vapor  calculation,  E,.(j)  and  Fu(j)  are  calculated  moving  down  the  grid, 

n  +  l  1 

and  6    (j)  is  calculated  moving  up  the  grid.   However,  8    (1)  must  be  known  in  advance. 

3. 6  Solution  of  the  Energy  Equation  for  the  Soil 

Assuming  constant  heat  capacity  and  constant  density  of  the  soil,  the  energy  equation 
for  the  soil  is 


3T    .       9T 

9t     9z   v  g  9z  ' 
g   S   g 


(72) 


Using  a  development  similar  to  that  used  in  the  water  vapor  calculation,  letting  £=  -=-, 
and  assuming  the  thermal  conductivity  of  the  soil  is  constant,  we  arrive  at  an  equation  of 
the  form 


"Ag(J)Tg+1(J+1)  =  BgU)Tg+1(J)  '  VJ^U-1)  =  Dg(j) 


(73) 
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where 


AM)  =  104  a(j)  K 


J  (j)  =  1  +  2.25  x  104  a(j)  K 
&  g 


C  (j)  =  1.25  x  104  a(j)  K 


Yj)  -  Tg«)  * lo4 «(»  S 


Tn(j+1)  -  2.25  Tn(j)  +  1.25  Tn(j-1) 


We  use  the  algorithm 


Tg+1(j)  =  Eg(j)  Tgn+1(j-l)  +  Fg(j)  (74) 


which  leads  to  the  expressions 


V35  =  yl 

and 


Eg(j3  =  B  (j)-A  (j)E  (j+l)  (7Sa) 


Dg(J)+Ag(j)Fg(j+l) 
Fg0)  =  B  (j)-A  (j)E  (j+l)  (75b) 

where  integral  values  of  j  correspond  to  points  on  the  z  grid. 

Applying  the  boundary  condition  on  T  at  the  lower  boundary  of  the  soil  layer 

dT 

— S.  =  0 

g 
we  get  the  relationship  T   (j  )  =  T   (j  -1),  where  j   is  the  grid  point  at  the  boundary. 

o        o         o        6  o 

Consequently,  from  equation  (74)  we  see  that 

E  (j  )  =  1.0  and  F  (j  )  =  0.0 
g   g  g  V 

Using  this  scheme,  E  (j)  and  F  (j)  are  calculated  moving  from  the  bottom  of  the  soil 
layer  up  to  the  earth's  surface.  Then  the  new  temperatures  are  calculated  using  equation 
(74)  while  moving  from  the  earth's  surface  down  to  the  bottom  of  the  soil  layer.  The  tem- 
perature at  the  earth's  surface  at  the  next  time  step  is  determined  by  an  energy  balance 
calculation  at  the  earth's  surface. 

3.7  The  Energy  Balance  at  the  Earth's  Surface 

Before  the  soil  temperature  and  potential  temperature  calculations  can  be  carried  out, 
the  temperature  of  the  soil  at  the  earth's  surface,  and  the  potential  temperature  of  the  air 
at  the  earth's  surface  must  be  known  for  the  next  time  step.  The  temperature  at  the  surface 
at  the  next  time  step  is  determined  by  an  energy  balance  calculation  at  the  earth's  surface. 
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Figure  4:   The  energy  balance  at  the  earth's  surface 

Consider  the  flow  of  energy  into  and  out  of  the  control  volume  indicated  by  broken 
lines  in  Figure  4.  The  control  volume  covers  an  area  of  one  square  centimeter  of  the  earth's 
surface  and  extends  from  z  (3/2)  to  z(3/2).   The  heat  flux  out  of  the  control  volume  due  to 
convection  is 

-'  S  "H  If 

The  heat  flux  out  of  the  control  volume  due  to  conduction  of  heat  into  the  soil  is 

9T 

-p   C   K  T-&- 
g   gr   g  3zg 

The  heat  flux  out  of  the  control  volume  due  to  evaporation  of  water  is 

-l»«vU 

where  L  is  the  latent  heat  of  vaporization.   The  energy  flux  out  of  the  control  volume  due  to 
emission  of  longwave  radiation  is 

E  B 
s  s 

where  e   is  the  emissivity  of  the  earth's  surface,  and  B   is  the  black  body  radiation  cor- 
responding to  the  surface  temperature. 

The  energy  flux  into  the  control  volume  due  to  absorption  of  longwave  radiation  is 

e  R  . 
s  net 

where  it  is  assumed  that  the  absorptivity  equals  the  emissivity.   The  energy  flux  into  the 
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control  volume  due  to  absorption  of  solar  radiation  is 


a  R  , 
s  solar 


where  a   is  the  absorptivity  of  solar  radiation  at  the  surface. 
In  the  model,  a   is  specified  as  a  function  of  time. 

The  rate  of  change  of  the  energy  stored  in  the  control  volume  is  the  sum  of  these  fluxes . 

3T 


p  C   z(3/2)  +  p   C   z  (3/2) 
P  g  gr  g 


3T     n      v ,  99       „    v 

W  =  °   Cp  KH  87  +  pg  Cgr  g  3i 


g 


+  L  p  K„  ^  +  e  R   .  -  e   B  +a  R  .     (76a) 
V  3z    s  net    s   s    s   solar   *•   ' 


We  substitute  expressions  for  z(3/2),  z  (3/2),  and  the  first  derivatives  into  the  equation 

1  g 

letting  £  =  -. 


IP  C   800  (/1.25  -1.0)  +  p   C   8(/I725".  1.0) 
P  g  gr 


Tn+1(l)-Tn(l) 


At 


l/2p  Cp  1^(3/2) 


+  l/2p   C    K 
g   gr   g 


3n+1(2)-en+1(i)+en(2)-en(i) 

800(ln  1.25)/1.25 
Tn+1(2)-Tn+1(l)+Tn(2)-Tn(l) 


8(ln  1.25)/1.25 


n+l,0.   n+1,.-.   n,0.   n,.. 
+  1/2  L  PKv(3/2)  3 HhZ W+<\    (2)-q  (1) 

800(ln  1.25)/1.25 


+  e(R    -  B  )  +  a  R  , 

s   net    s     s  solar 


(76b) 


We  make  the  following  substitutions: 


3n+1(2)  =  EH(2)  6n+1(l)  +  FH(2) 


)n+1(l)  =  Tn+1(l) 


P(l) 


Tn+1(2)  =  E  (2)  Tn+1(l)  +  F  (2) 


Tn+1(l)  =  Tn+1(l) 


We  assume  that  qn+  (2)  =  q"(2)  and  q"+  (1)  =  q°(l)  since  q  varies  slowly.   Collecting  terms 
with  T   (1)  on  the  left  hand  side  of  the  equation,  we  get 


„n+l 


At 


T        (1){(100  p   C     +   p      C      )8(/iT25~  -   1.0) 

P  g     gr  1600(ln/T725)    1.25 
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P     K 
O 


(EH(2)  -  1.0)  (^y)   P  Cp  K^(3/2)  +  100  Pg  Cgr  Kg  (Eg(2)  -  1.0) 


} 


Tn(l)  8(100  p  C   +  p   C   )(/l.25  -  1.0)  +  — ■ 

P    §  gr  1600  (In  1.25)/1.25 


{pC   K^(3/2) 


FH(2)  +  en(2)  -  en(i) 


+  100  p   C   K 
g   gr   g 


+  2  L  p  Ky(3/2) 


F  (2)  +  Tn(2)  -  Tn(l) 


qn(2)  -  qn(l) 


Let 


+  At|e(R    -B)+a   R   , 
1  s  net    s'    s  solar 


(76c) 


and 


Also  let 


'si 


's2 


1600  (In  1.25)/1.25 

At 


800(p(l)  C  +  0.01  p  C  X^l-25  -  1.0) 
P        g  gr 


D  =  1  -  C  ,  C  „ 

si   s2 


P  K 
o 


p(3/2)  C   K^(3/2)(EH(2)  -  1.0)  (^y) 


+  100  p   C   K   (E  (2)  -  1.0) 
g  gr  g   g 


n+1 


Solving  for  T   (1),  we  get 


Tn+1(l) 


T  (1)  +  C 


s2 


+100  p   C   K 
g   gr  J 


q"(2)  -  qn(l) 


Csl{p(3/2)  Cp  K^(3/2) 


FH(2)  +  6n(2)  -  9n(l) 


F  (2)  +  Tn(2)  -  Tn(l)l  +  2  L  p(3/2)  1^(3/2) 

-  J  n 


„n+l 


} 


e  (R    -  B  )  +  a  R  . 
s^  net    s     s  sol 


ar 


After  T   (1)  has  been  calculated,  then 

K 


/D 


and 


)n+1(i)  =Tn+1(i)  (p2_) 


Tn+1(l)  =  Tn+1(l) 


(77) 


3.8  The  Specific  Humidity  at  the  Earth's  Surface 

When  the  relative  humidity  is  specified  as  a  constant  at  the  earth's  surface,  the 
specific  humidity  is  calculated  for  the  next  time  step  using  the  surface  temperature  result- 
ing from  the  energy  balance  calculation.   When  the  latent  heat  flux  at  the  earth's  surface 
F  is  specified  as  a  function  of  time,  the  specific  humidity  is  calculated  using  the 
relationship 


Fq  =  "P  L  "V 


IS. 
9z 


(78a) 
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Expressing  this  in  finite  difference  form,  we  get 

n+1 


n+1 


F   =  -p(l/2)  L  ICM/2)  -^ (2)  -  q    (1)_ 

q  800(ln  1.25)/1.25 


(78b) 


Making  the  substitution  q    (2)  =  E  (2)  q   (1)  +  F  (2),  and  solving  for  q   (1),  we  get 


F  (2)  +  800(ln  1.25)/T25  F  /[L  p(l/2)  K  (1/2)1 
qu  A(l)  =  — H_Jz » -1 


(79) 


1  -  Ev(2) 


The  decision  regarding  which  boundary  condition  to  use  depends  upon  the  surface  conditions 
and  upon  the  available  data. 

3 .9  The  Test  for  Saturation 

Let  T(j)  and  q(j)  be  the  tentative  values  of  temperature  and  specific  humidity  of  water 
vapor  resulting  from  the  calculations  described  above,  and  let  S(T)  be  the  saturation  density 
of  water  vapor  at  the  temperature  T.   It  is  possible  to  determine  whether  or  not  the  air  is 
saturated  by  comparing  p(j)q(j)  with  S(T(j)).   If  p(j)q(j)  <  S(T(j)),  the  air  is  unsaturated, 
and  no  adjustment  of  values  is  necessary. 

If  p(j)q(j)  >  S(T(j)),  the  air  is  supersaturated,  and  the  values  of  q  and  T  must  be  ad- 
justed to  return  the  system  to  saturation  conditions. 


pq 
pq~ 


•q 


+  l 


1 

Ad) 

3, 

2 

hm 

T  T 


Figure  5:  The  adjustment  for  supersaturation 


3S 
Let  S'(T)  =  — ,  and  we  approximate  S'(T)  at  point  3  in  Figure  5  by 

a  I 

.n+1 


s,^(j))=S(T-(j))  -S(1(J)) 
Tn  '(j)  -  T(j) 


(80) 


From  an  energy  balance  consideration,  we  get 


or 


C   AT  =  -L  Aq 
P  H 

Tn+1(j)  -TO)  =  -i-(qn+1(j)  -5(3)) 


(81) 
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Saturation  conditions  exist  at  T 


n+1 


hence 


S(T    (j))  =  P  q    (j) 

We  substitute  these  relationships  into  equation  (80)  and  solve  for  q   (j),  which  leads  to 
the  result 

qn+l(j]  =  sftCj))  +  (L/Cp)  S'(^(3))q(3)  (82) 

P(j)  +  (L/C  j  S'(T(j)) 

Now  that  q   (j)  is  known,  T   (j)  can  be  calculated  using  equation  (81).   In  the  nu- 
merical model,  values  of  S  and  S'  are  tabulated  for  temperature  intervals  of  0.5°K  in  the 
temperature  range  of  223°K  to  320°K.   When  T  lies  between  temperatures  used  in  the  tabula- 
tion, values  of  S  and  S'  at  the  next  higher  tabulated  temperature  are  used.   The  water  vapor 
is  assumed  to  condense  in  the  form  of  fog  or  stratus,  and  the  liquid  water  is  allowed  to 
remain  in  the  atmosphere  until  it  has  evaporated. 

3.10  Solution  of  the  Longwave  and  Solar  Radiation  Equations 

For  the  layer  between  two  grid  points,  assume  a  linear  relationship  between  B  and  id. 
For  the  layer  between  grid  points  j  and  j+1,  we  have 


and 


B.    .    -   B. 

1.      +    -J*i 1    (CD     -     <D.) 

J  Mj+1     "    0).  J 


(83) 


dB   = 


-  J+1 


J 


(D  .      ,      -ID. 

3+1  J 


d(D 


(84) 


Integration  can  now  be  carried  out  analytically  between  adjacent  grid  points  for  the  long- 
wave radiation  calculation.   For  example: 


(D.   , 


t(id)  dB(co) 


ID.  .  B.  . 
J+1       J+1 


(D.  .. 


*L  z       a,    exp(-B,(D)dt 


j  k=l 


or 


j+1T(U)dB(U)  =  -   J*1  "   j  I 


\    r 


(D.   .   -  (D.  . 

j+1    j   k=l 


exp(-6k  (d   )-exp(-6k  u.)| 

J 


The  integrals  in  equation  (31)  can  now  be  evaluated  by  integrating  over  the  layers  be- 
tween grid  points  and  summing  over  the  layers  included  in  the  limits  of  the  integral.   Let 
j   be  the  grid  point  corresponding  to  the  reference  height,  and  let  j    be  the  grid  point 
at  the  top  of  the  grid  where  <d  =  0.   At  the  earth's  surface,  j  =  1  where  id  =  W.   After  the 
integrals  have  been  replaced  by  summations,  we  get 

9R 


net 


3(D 


=  e   B 
s   s 


6 

I 

k=l 


ak  Bk  exp 


-6k(o>1  -  u.  ) 
r  _ 


-  B.   I   a      B,  exp 
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The  absorptivity  of  water  vapor  is  affected  by  changes  in  pressure.   This  effect  is  in- 
cluded in  the  calculation  of  the  optical  thickness.   An  effective  absorbing  mass  is  defined 
as 


pq(P/P0) 
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When  this  is  included  in  the  definition  of  optical  thickness,  we  get 
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We  approximate  the  integral  by  calculating  the  change  in  optical  thickness  Ato .  between  grid 


points  z.  and  z.  ,,  then  summing. 
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where  pq(P/P  )     is  the  average  value  between  z.  and  z.  ..   From  the  definition  of  the 
o  6  j      3+1 

optical  thickness  we  see  that 
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When  this  is  substituted  into  equation  (29a), we  see  that  the  radiative  flux  divergence  at 
the  grid  point  j  is 


3R   .                n,..  1/2   3R   . 
,   net                 rP(jJ\     r      neti 
(^ )■    =  -P(J)  q(j)  (^-)  (^ )■ 


(89) 


29 


When  equations  (38b)  and  (88)  are  substituted  into  equation  (37),  we  see  that  the  divergence 
of  solar  radiation  at  the  j   grid  point  is 
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3.11  General  Plan  of  Calculation 

The  calculation  begins  with  an  assumed  or  a  known  set  of  initial  conditions.   Initially, 
the  absolute  temperature,  specific  humidity,  and  components  of  the  wind  must  be  specified  at 
each  grid  point  in  the  atmospheric  region,  and  the  absolute  temperature  must  be  specified  in 
the  soil  region.   The  horizontal  gradient  of  temperature  must  also  be  specified  at  each  grid 
point  in  the  atmospheric  region  when  advection  is  included  in  the  calculation.   The  potential 
temperature  and  atmospheric  pressure  are  calculated  at  grid  points  given  the  atmospheric 
pressure  at  the  earth's  surface  and  the  temperature  and  specific  humidity  data. 

The  calculation  moves  forward  in  time  by  first  calculating  a  new  temperature  at  the 
earth's  surface  using  an  energy  balance  calculation.   The  new  values  for  the  potential 
temperature  in  the  atmosphere  and  the  absolute  temperature  in  the  soil  are  then  calculated 
for  the  next  time  step  using  the  new  surface  temperature.   Several  methods  may  be  used  to 
determine  a  new  value  for  the  specific  humidity  at  the  earth's  surface  depending  upon  the 
boundary  conditions.   Once  the  specific  humidity  is  determined  at  the  surface  for  the  next 
time  step  using  one  of  the  methods  described  in  section  2.2,  new  values  are  calculated  for 
the  specific  humidity  at  grid  points  in  the  atmospheric  region. 

After  the  new  temperature  and  specific  humidity  values  are  calculated,  a  test  is  made 
to  determine  whether  the  atmosphere  is  saturated.  If  it  is  supersaturated,  adjustments  are 
made  in  the  temperature  and  specific  humidity  to  bring  the  system  to  saturation  conditions 
accompanied  by  the  condensation  of  water  vapor  in  the  form  of  fog  or  stratus.  Conversely, 
if  the  atmosphere  is  unsaturated  and  liquid  water  is  present  in  the  form  of  fog  or  stratus, 
evaporation  is  allowed  to  occur  with  changes  in  temperature  and  specific  humidity  so  as  to 
maintain  saturation.  With  these  new  values  for  the  temperature  and  specific  humidity,  the 
process  of  computing  values  for  the  next  time  step  is  repeated. 

Various  components  of  the  energy  budget  are  calculated  to  provide  a  more  detailed  anal- 
ysis of  the  processes  involved  in  the  calculation.  The  components  of  the  energy  balance 
calculation  at  the  earth's  surface  are  also  computed  separately.  Although  the  prediction  of 
accurate  temperature  and  specific  humidity  values  is  the  primary  goal,  the  accurate  predic- 
tion of  components  of  the  energy  budget  must  also  be  our  goal.   Since  the  numerical  model  is 
the  result  of  combining  models  of  the  various  energy  transfer  processes  in  the  atmosphere, 
the  total  model  will  be  inaccurate  if  any  one  of  these  processes  is  modeled  inaccurately. 
It  is  also  possible  that  the  errors  in  the  models  of  the  energy  transfer  processes  may  be 
compensating,  resulting  in  a  false  sense  of  accuracy  in  the  overall  calculation.   Consequently 
the  accuracy  of  numerical  calculation  should  be  analyzed  in  terms  of  both  the  overall  ac- 
curacy and  the  accuracy  of  various  processes  involved  in  the  calculation. 

The  terrestrial  radiation  calculation  is  performed  every  fifteen  minutes  rather  than  at 
every  time  step.   This  is  done  because  the  radiation  flux  and  flux  divergence  values  vary 
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slowly,  thus  permitting  a  considerable  savings  in  computing  time  by  not  having  to  perform 
the  calculation  at  every  time  step.   The  solar  radiation  calculation  is  performed  at  every 
time  step,  however. 

The  horizontal  components  of  the  wind  field  are  specified  for  the  duration  of  the  cal- 
culation.  The  horizontal  gradient  of  temperature  in  the  direction  of  the  wind  is  also 
specified  for  the  duration  of  the  calculation  when  advection  is  included.   New  wind  data 
and  horizontal  temperature  gradient  data  are  read  in  at  regular  intervals.   The  wind  data 
are  used  in  the  calculation  of  the  eddy  diffusion  coefficient  for  heat  and  in  the  calculation 
of  the  advection  term  in  the  energy  equation. 

In  the  preparation  of  initial  conditions  and  wind  data  it  was  necessary  to  interpolate 
between  values  of  the  observational  data  available.   A  spline  fit  was  used  whenever  inter- 
polation was  necessary.   According  to  this  method,  a  smooth  curve  with  continuous  slope  and 
curvature  is  generated  which  passes  through  the  given  data  points.   Interpolation  is  done  by 
merely  reading  points  off  the  curve.   The  spline  fit  is  accomplished  by  connecting  each 
pair  of  adjacent  points  with  a  section  of  a  third  degree  polynomial,  matching  up  the  sections 
so  that  the  first  and  second  derivatives  are  continuous  at  each  point.   At  the  end  points 
an  additional  requirement  is  made  that  the  second  derivative  be  a  linear  extrapolation  of 
the  value  at  the  two  adjacent  points. 

4.0   RESULTS 

4 . 1   Introduction 

In  the  area  of  numerical  modeling,  the  presentation  of  results  is  complicated  by  the 
volume  of  data  output.   In  an  effort  to  be  complete  yet  still  maintain  clarity,  results  are 
presented  in  graphical  and  tabular  form.   Results  are  presented  in  terms  of  the  overall  cal- 
culation first,  then  the  various  components  of  the  calculation  are  analyzed  separately. 

By  comparing  computed  values  of  the  temperature  and  moisture  fields  and  the  components 
of  the  energy  budget  with  the  available  data,  the  accuracy  of  the  model  can  be  checked.   This 
is  an  important  step,  for  the  limitations  of  the  model  must  be  known  if  the  results  are  to 
be  trusted  for  accuracy.   Ideally,  the  model  should  be  tested  under  a  variety  of  conditions 
to  test  its  versatility.   We  are  limited  in  our  choice  of  conditions  by  the  data  which  are 
available  for  comparison,  nevertheless  it  would  be  nice  if  the  model  could  be  tested  under 
conditions  of  strong  winds  and  light  winds  and  under  conditions  where  advection  is  not  sig- 
nificant and  under  conditions  where  it  is  significant.   The  O'Neill,  Nebraska  data  (Lettau 
and  Davidson,  1953)  and  the  Davis,  California  data  (Pruitt,  Morgan,  and  Lourence,  1968) 
provide  a  variety  of  conditions  for  comparison.   Estoque's  calculations  (1963)  are  also  a- 
vailable  for  comparison  with  the  O'Neill  data,  so  the  model  can  also  be  compared  with  another 
numerical  model,  which  is  of  a  very  specialized  nature. 

The  O'Neill  test  site  is  located  on  flat  terrain  with  a  nearly  uniform  grass  cover  and 
is  devoid  of  trees.   Measurements  were  made  under  strong  southerly  wind  conditions.   The 
primary  obstruction  to  general  southerly  flow  is  found  in  the  regions  to  the  west,  south, 
and  east  of  the  site  and  are  known  as  the  Sand  Hill  and  Loess  Hill  country.   These  regions 
are  rolling  hills  which  usually  are  not  more  than  100  ft.  high.   In  general  the  rolling  hills 
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become  less  pronounced  as  the  test  site  is  approached,  so  the  area  within  a  10  to  20  mile 
radius  of  the  site  is  almost  free  of  any  hilly  obstructions.   Advection  is  considered  a 
small  effect  in  this  region,  which  was  the  primary  reason  for  its  choice  as  a  test  site. 

The  Davis  site  is  located  on  flat  terrain  about  15  miles  from  the  coast  range,  45  miles 
from  the  Sierras,  and  17  miles  west  south  west  of  Sacramento  in  the  Central  Valley  of  Cal- 
ifornia.  The  site  is  in  a  general  agricultural  area  removed  from  any  significant  buildings 
or  changes  in  landscape.   Prevailing  winds  were  from  the  south  west  or  north  and  were  mod- 
erate to  light  during  the  time  measurements  were  made.   Advection  has  been  shown  to  be  sig- 
nificant in  this  region  by  Fosberg  and  Schroeder  (1966). 

The  calculations  under  consideration  are: 

(1)  O'Neill,  Nebraska;  1230  August  24,  1953  to  1230  August  25,  1953.  (LST) 

(2)  Davis,  California;  0100  May  4,  1967  to  0100  May  5,  1967. (LST) 

These  calculations  are  henceforth  referred  to  as  the  O'Neill  calculation  and  Davis  cal- 
culation, respectively.   These  calculations  are  conducted  first  neglecting  the  effects  of  ad- 
vection.  The  advection  term  is  then  included  in  the  energy  equation,  and  the  magnitude  of 
the  horizontal  temperature  gradients  necessary  to  make  the  calculated  temperature  values 
agree  with  observed  values  are  estimated. 


4.2  Results  of  the  O'Neill  Calculation 

The  following  data  are  used  in  the  O'Neill  calculation: 
AT   =  30.0  sec  Latitude  =  42.470  degrees  N 

C    =  0.360  cal/gm  deg  Solar  Declination  =  11.030  deg 

K    =  0.0015  cm2/sec  Sunrise  =  0551  LST 

^3  2 

p    =1.40  gm/cm  I  =  1.94  cal/cm  min 

KH*  =  6.00  x  105  cm2/sec  e   =  0.90 

C:   =  50.0  U  =  0.861 

C2   =  6.00  x  106  cm2/sec  d  =  8.00  cm 

Surface  Pressure  =  0.940  atm 

The  variable  y      is  set  equal  to  0.65  °C/km  above  200  meters  height,  and  it  decreases  to  a 

value  of  zero  at  100  meters.   The  specific  humidity  is  specified  at  the  earth's  surface  as 

a  function  of  time  in  the  calculation. 

Temperature 

The  computed  temperature  field  for  the  24-hour  period  beginning  at  1230  CST  on  August 
24,  1953  is  shown  in  Figures  6  and  7.   Comparison  between  observed  and  predicted  temperature 
profiles  is  shown  in  Figures  8  through  12.   The  general  features  of  the  predicted  temperature 
variations  agree  well  with  the  observed  variations,  however  there  is  some  discrepancy  be- 
tween the  two  at  2130  CST.   The  numerical  model  predicts  a  slight  warming  trend  in  the  region 
between  20  and  100  meters  height  from  2130  to  2230,  but  this  trend  was  not  present  in  the 
observed  temperature  field.   This  appears  as  a  large  "bump"  in  the  302°K  contour  line  at 
2100  CST  in  Figure  7. 

The  cause  of  the  warming  trend  is  apparent  when  we  look  at  the  variations  in  the 
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Richardson  number  and  eddy  diffusion  coefficient  for  heat  in  this  region.   The  Richardson 
number  reaches  a  local  maximum  at  a  height  of  26  meters  of  0.269  at  2130  CST  compared  with 
values  of  0.067  and  0.073  a  few  meters  above  and  below  that  level.   The  result  is  that  the 
eddy  diffusion  coefficient  has  a  local  minimum  at  a  height  of  20  meters,  thus  restricting 
the  downward  flow  of  heat  by  convection  out  of  the  region. 

The  Richardson  number  has  a  local  minimum  at  a  height  of  137  meters  with  a  correspond- 
ing local  maximum  in  the  eddy  diffusion  coefficient  for  heat  at  the  same  height.   Looking  at 
the  potential  temperature  profile  at  2130  CST  in  Figure  13,  we  see  that  the  potential  tem- 
perature increases  with  height.   Consequently,  energy  is  moved  downward  by  convection  into  a 
region  where  there  is  a  restriction  in  the  convective  heat  flux,  resulting  in  an  accumulation 
of  energy  and  a  warming  trend  in  that  region.   This  would  not  have  occurred  if  the  Richardson 
number  had  not  varied  so  greatly  over  such  a  small  range  of  heights.   This  problem  may  be 
overcome  by  smoothing  out  the  Richardson  number  by  an  averaging  process  in  the  first  100 
meters  or  so  of  the  atmosphere. 

Temperature  profiles  are  plotted  for  6-hour  intervals  beginning  at  1230  CST  in  Figures 
8  through  12.   The  observed  temperature  profiles  and  those  predicted  by  Estoque  are  also 
plotted.   Figure  8  shows  the  initial  conditions  for  the  calculation.   In  Figures  9  and  10  we 
see  that  the  predicted  temperatures  are  slightly  higher  than  those  observed,  but  the  pre- 
dicted profiles  are  of  the  same  general  shape  as  those  observed.   This  is  not  always  the 
case  with  Estoque' s  temperature  profiles. 

Figure  11  shows  very  good  agreement  between  the  predicted  and  observed  temperatures  up 
to  a  height  of  150  meters.   The  calculation  predicts  a  temperature  inversion  which  was  ob- 
served, although  it  occurred  at  a  different  height  than  predicted.   Estoque 's  calculation 
shows  a  noticeable  variation  from  the  other  two  profiles. 

Figure  12  shows  excellent  agreement  between  the  observed  and  predicted  temperature 
profiles.   Estoque 's  calculations  again  are  not  as  accurate  as  those  presented  here. 

The  Specific  Humidity 

The  computed  moisture  field  for  the  24-hour  period  is  shown  in  Figures  15  and  16,  and 
comparison  between  predicted  and  observed  moisture  profiles  is  shown  in  Figures  17  through 
21.  The  general  features  of  the  moisture  calculation  agree  fairly  well  with  observed  values, 
but  the  results  are  not  as  good  as  those  for  the  temperature  calculation.  As  in  the  tem- 
perature calculation,  there  is  a  discrepancy  in  the  moisture  calculation  at  2130  CST  which 
appears  as  a  peak  in  the  11.25  gm/kgm  contour  line  in  Figure  16.   The  specific  humidity  be- 
tween 20  meters  and  100  meters  height  decreases  slightly  between  2130  and  2230  rather  than 
continuing  to  increase  as  was  observed.   The  cause  is  the  same  as  that  found  in  the  temper- 
ature calculation.   The  assumption  was  made  that  the  eddy  diffusion  coefficient  for  water 
vapor  equaled  the  eddy  diffusion  coefficient  for  heat.   Consequently,  when  the  eddy  diffusion 
coefficient  for  heat  had  a  local  minimum  at  20  meters,  the  eddy  diffusion  coefficient  for 
water  vapor  also  had  a  local  minimum.   The  result  is  that  the  upward  flux  of  water  vapor  is 
restricted  at  a  height  of  20  meters,  but  it  is  not  restricted  at  higher  levels.   Consequently, 
water  vapor  moves  upward  out  of  the  region  above  20  meters  with  the  result  that  the  specific 
humidity  decreases  in  this  region. 
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The  specific  humidity  profiles  are  plotted  for  6-hour  intervals  in  Figures  17  through  21. 
The  observed  profiles  and  those  predicted  by  Estoque  are  also  plotted.   There  is  good  agree- 
ment between  predicted  and  observed  moisture  profiles  for  the  first  12  hours  of  the  calcul- 
ation, and  the  predicted  profiles  show  the  same  general  variations  as  the  observed  profiles 
for  the  remainder  of  the  calculations.   A  sharp  decline  in  specific  humidity  is  predicted 
in  Figure  21  by  the  numerical  model,  and  a  sharp  decline  was  observed,  although  the  height 
at  which  the  decline  began  occurred  at  a  slightly  different  height. 

There  are  several  reasons  for  the  discrepancies  between  predicted  and  observed  moisture 
profiles.   The  simplified  boundary  conditions  for  the  water  vapor  calculation  at  the  earth's 
surface  is  the  primary  reason  for  discrepancies  near  the  earth's  surface.   The  observed 
moisture  profiles  showed  large  variations  in  the  specific  humidity  above  1000  meters.   For 
example,  at  a  height  of  2000  meters  the  specific  humidity  varied  between  10.4  gm/km  and  2.5 
gm/km  during  the  24-hour  period  observed.   At  a  height  of  1000  meters,  the  specific  humidity 
varied  between  11.6  gm/kgm  and  9.5  gm/kgm  in  two  hours.   Variations  such  as  these  are  pro- 
bably attributable  to  advection,  which  was  not  included  in  the  water  vapor  calculation. 

The  assumption  was  made  that  the  eddy  diffusion  coefficients  for  heat  and  water  vapor 
were  equal.   Experimental  evidence  supports  this  assumption  near  the  earth's  surface  (Craw- 
ford, 1965a), but  this  assumption  may  be  invalid  at  higher  levels. 

The  Energy  Balance  at  the  Earth's  Surface 

The  predicted  energy  budget  constituents  at  the  earth's  surface  for  the  O'Neill  cal- 
culation are  shown  in  Table  3  along  with  the  average  of  the  values  calculated  by  Suomi  and 
Lettau  (Lettau  and  Davidson,  1957).   The  predicted  sensible  heat  flux  F„  agrees  well  with 
Suomi  and  Lettau' s  values,  although  the  predicted  values  are  somewhat  larger  between  2030 
and  0630  CST.   The  good  agreement  between  these  values  and  the  good  agreement  between  the 
predicted  and  observed  temperature  profiles  indicate  that  the  model  for  the  eddy  diffusion 
coefficient  is  quite  good.   The  eddy  diffusion  coefficient  would  be  expected  to  increase 
with  height  until  it  reached  a  maximum  at  a  height  of  about  a  hundred  meters,  after  which 
it  would  be  expected  to  decrease  gradually.   The  eddy  diffusion  coefficient  corresponding  to 
1630  on  August  24  shows  this  variation.   The  eddy  diffusion  coefficient  increases  rapidly 

f\  9        1 

with  height,  reaching  a  maximum  value  of  2.67  x  10  cm  sec   at  a  height  of  67  meters  and 

5   2    - 1 
decreases  gradually  to  a  value  of  1.2  x  10  cm  sec   at  2100  meters.   The  large  values  for 

the  eddy  diffusion  coefficient  reflect  the  increased  turbulent  mixing  in  the  atmosphere  re- 
sulting from  strong  winds. 

Comparing  the  predicted  latent  heat  flux  F  with  Suomi  and  Lettau' s  values,  we  see  that 
the  agreement  is  not  very  good.   Positive  values  as  large  as  0.107  ly  min   are  predicted  be- 
tween 2030  and  1630  while  the  values  calculated  by  Suomi  and  Lettau  are  negative  during  this 
period.   Ilalstead  also  calculated  the  latent  heat  flux  at  O'Neill  (Lettau  and  Davidson,  1957) 
beginning  at  0030  on  August  24,  1953.   His  values  show  better  agreement  with  the  predicted 
values.   His  values  for  the  latent  heat  flux  were  also  positive,  and  they  ranged  from  0.068 
ly  min   at  0030  to  0.285  ly  min"  at  0830.   The  large  discrepancies  shown  here  between  the 
various  calculations  make  it  difficult  to  determine  the  accuracy  of  the  latent  heat  flux 
calculation.   Considering  the  accuracy  of  the  predicted  temperature  and  moisture  profiles, 
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Figure  6.    The  computed  temperature  field  up  to  0.50  bar  August  24,  1953  at  O'Neill,  Nebraska. 
First  solid  line  is  280°  K  with  2°  intervals. 


2000.0  rr~i — i,.i..j.-i-^t-^-^->-^ — >—  -^ i i--L — L-J- *— •  i — L i — < — i — 


1000.0-- 


500.0 


J 1 — K— I — H— I — I — H^H 


1230      1430      1630      1830      2030      2230      0030      0230     0430      0630      0830      1030       1230 
Figure  7.    The  computed  temperature  field  up  to  2000  meters  for  August  24,  1953  at  O'Neill,  Nebraska. 
First  solid  line  is  290°K  with  a  2°  interval. 
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Figure  8.    Temperature  profiles  1230  CST  on  August  24,  1953  at  O'Neill,  Nebraska. 
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Figure  9.    Temperature  profiles  1830  CST  on  August  24,  1953  at  O'Neill,  Nebraska. 
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Figure  10.    Temperature  profiles  0030  CST  on  August  25,  1953  at  O'Neill,  Nebraska. 
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Figure  11.    Temperature  profiles  0630  CST  on  August  25,  1953  at  O'Neill,  Nebraska. 
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Figure  12.    Temperature  profiles  1230  CST  on  August  25,  1953    at  O'Neill,  Nebraska. 
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Figure  13.    The  computed  potential  temperature  field  up  to  0.50  bar  for  August  24,  1953  at  O'Neill, 


Nebraska. 
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Figure  14.    The  computed  potential  temperature  field  up  to  2000  meters  for  August  24,  1953  at 
O'Neill,  Nebraska. 
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Figure  15.    The  computed  specific  humidity  (gm/kgm)  field  up  to  0.50  bar  for  August  24,  1953  at 

O'Neill,  Nebraska. 
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Figure  16.    The  computed  specific  humidity  (gm/kgm)  field  up  to  2000  meters  for  August  24,  1953 
at  O'Neill,  Nebraska. 
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Figure  17.    Specific  humidity  profiles  1230  CST  on  August  24,  1953  at  O'Neill,  Nebraska. 
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Figure  18.    Specific  humidity  profiles  1830  CST  on  August  24,  1953  at  O'Neill,  Nebraska. 
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Figure  19.    Specific  humidity  profiles  0030  CST  on  August  25,  1953  at  O'Neill,  Nebraska. 
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Figure  21.    Specific  humidity  profiles  1230  CST  on  August  25,  1953  at  O'Neill,  Nebraska. 
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Figure  22.    The  computed  relotive  humidity  field  up  to  0.50  bar  for  August  24,  1953  at  O'Neill, 
Nebraska. 
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Figure  23.    The  computed  relotive  humidity  field  up  to  2000  meters  for  August  24,  1953  at  O'Neill, 
Nebraska. 
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it  is  not  likely  that  the  errors  are  very  large  in  the  latent  heat  flux  calculation. 

The  predicted  heat  flux  into  the  soil  F   agrees  well  with  Suomi  and  Lettau's  values. 

gr 

Agreement  is  particularly  good  between  2030  and  0830  CST.   The  thermal  conductivity  of  the 

soil  was  assumed  constant  throughout  the  soil  region,  and  this  assumption  may  account  for 

some  of  the  discrepancy  between  the  two  sets  of  values. 

Comparison  of  the  two  sets  of  radiation  flux  data  also  shows  very  good  agreement.   This 

is  particularly  encouraging  because  the  net  values  for  the  radiation  flux  F   ,  are  the  re- 

rad 

suit  of  two  radiation  calculations—the  terrestrial  radiation  calculation  and  the  solar  ra- 
diation calculation—and  errors  in  either  one  of  the  calculations  would  result  in  errors  in 
the  net  value. 

4.3  Results  of  the  Davis  Calculation 

The  following  data  are  used  in  the  Davis  calculation: 
AT   =30.0  sec  Latitude  =  38.540  degrees  N 

C    =  0.525  cal/gm  deg  Solar  Declination  =  15.90  degrees 

K    =  0.0029  cm2/sec  Sunrise  =  0511  LST 

3  2 

p    =1.30  gm/cm  I  =1.94  cal/cm  min 

K  *   =  6.00  x  105  cm2/sec  e   =  0.90 

n  S 

C    =  50.0  y  =  0.804 

&  9 

C2   =  6.00  x  10   cm  /sec  d  =  10.6  cm 

Surface  Pressure  =  1.012  atm 

Temperature 

The  predicted  temperature  field  for  the  24-hour  period  beginning  at  0100  PST  on  May  4, 
1967  is  shown  in  Figures  24  and  25.  Comparison  between  observed  and  predicted  temperature 
profiles  is  shown  in  Figures  26  through  30.  The  predicted  temperature  field  shows  the  same 
general  variations  as  the  observed  temperature  field. 

When  we  refer  to  the  temperature  profiles  plotted  for  6-hour  intervals  in  Figures  26 
through  30,  we  see  that  there  is  quite  a  variation  in  agreement  between  predicted  and  ob- 
served profiles.  There  is  very  good  agreement  between  the  two  temperature  profiles  in 
Figure  28,  however  there  are  discrepancies  as  large  as  6°C  in  Figures  27  and  29.  There  are 
two  primary  reasons  for  these  discrepancies.  Advection  is  known  to  be  a  contributing  factor 
in  causing  some  of  the  discrepancy,  and  this  probably  accounts  for  most  of  the  discrepancy 
in  Figure  29.   Cooler  air  usually  moves  into  the  Davis  area  in  the  afternoon  from  the  San 
Francisco  Bay  Area,  which  lies  about  60  miles  southwest  of  the  Davis  site.   This  would  ac- 
count for  the  cooling  observed  in  the  actual  temperature  field  between  1300  and  1900  PST. 
Nevertheless,  one  would  still  expect  a  larger  drop  in  the  predicted  temperature  between 
1900  PST  and  0100  PST  than  that  shown.  The  problem  seems  to  be  that  too  little  energy  is 
being  transferred  from  the  atmosphere  to  the  soil  under  conditions  of  strong  thermodynamic 
stability  near  the  earth's  surface.  Referring  to  Table  4  we  see  that  very  little  energy  is 
transferred  to  the  soil  by  convection  during  the  times  when  the  surface  is  cooling.   Since 
the  potential  temperature  increases  vertically,  heat  is  moved  downward  into  the  region  near 
the  earth's  surface,  resulting  in  a  prediction  of  higher  temperatures  than  those  observed. 
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The  problem  is  related  to  the  calculation  of  the  eddy  diffusion  coefficient.   Crawford 
(1965)  noted  that  accurate  calculation  of  the  convective  heat  flux  near  the  earth's  surface 
is  particularly  difficult  under  conditions  of  strong  thermodynamic  stability.   The  value  of 
<)>(R.)  is  uncertain  for  Richardson  numbers  greater  than  0.1  because  of  the  widely  scattered 
data  points  used  in  the  determination  of  <j>(R.)-   Using  the  expression  for  <|>(R.)  described 
above  for  the  forced  convection  regime  results  in  negligible  convective  heat  fluxes  for 
Richardson  numbers  greater  than  0.3.   In  the  Davis  calculation  the  computed  Richardson 
number  at  one  meter  reached  values  greater  than  1.0  between  1900  and  0100  PST,  which  accounts 
for  the  small  values  for  the  predicted  convective  heat  flux.   The  Richardson  number  depends 
upon  the  gradient  of  the  wind  speed.   Accurate  wind  data  are  important  for  the  first  few 
meters  of  the  atmosphere  especially  under  conditions  of  light  wind.   Uncertainties  in  the 
values  of  the  wind  data  may  account  for  the  large  Richardson  numbers  near  the  earth's  surface 
in  the  calculation. 

Specific  Humidity 

The  predicted  moisture  field  is  shown  in  Figures  33  and  34  along  with  the  moisture 
profiles  which  are  plotted  for  6-hour  intervals  in  Figures  35  through  39.   The  predicted 
moisture  field  shows  the  general  variations  observed  in  the  actual  moisture  field.   Agree- 
ment between  the  observed  and  predicted  moisture  profiles  is  very  good  below  250  meters. 
The  noticeable  variations  in  the  observed  moisture  field  above  this  level  can  probably  be 
attributed  to  advection  effects. 

The  Energy  Balance  at  the  Earth's  Surface 

The  predicted  energy  budget  constituents  at  the  earth's  surface  for  the  Davis  calcul- 
ation are  shown  in  Table  4  along  with  the  values  measured  by  Pruitt,  Morgan,  and  Lourence 
(1968).   The  predicted  sensible  heat  flux  agrees  very  well  with  the  measured  values,  although 
the  magnitudes  of  the  predicted  values  are  somewhat  smaller  than  observed  values  during  the 
period  from  1700  to  0100  PST. 

The  eddy  diffusion  coefficient  at  1700  on  May  9,  1967  reaches  a  maximum  value  of  4.66 

5   2-1  4   2-1 

x  10  cm  sec   at  a  height  of  67  meters  and  decreases  to  a  value  of  3.73  x  10  cm  sec   at 

a  height  of  2100  meters.   The  eddy  diffusion  coefficient , therefore,  behaves  very  similarly 

to  that  in  the  O'Neill  calculation.   The  smaller  values  for  the  eddy  diffusion  coefficient 

for  the  Davis  calculation  indicate  the  decrease  in  turbulent  mixing  resulting  from  lighter 

winds . 

The  latent  heat  flux  at  the  earth's  surface  was  specified  for  this  calculation,  so  agree- 
ment between  the  two  sets  of  values  is  assured.  The  heat  flux  into  the  soil  shows  fair  agree- 
ment with  the  measured  values.  The  agreement  was  not  as  good  in  this  case  as  it  was  in  the 
O'Neill  calculation. 

The  predicted  net  radiation  flux  agrees  very  well  with  the  measured  values,  as  was  also 
the  case  in  the  O'Neill  calculation.  This  demonstrates  further  the  accuracy  of  the  ter- 
restrial and  solar  radiation  calculations. 


46 


rt 

<u 

0 

X, 

+j 

+-> 

03 

ri 

•H 

e 

w 

U 

+j 

o 

c 

<4-l 

(U 

•H 

3 

i — 1 

+-> 

ol 

•H 

U 

■M 

(/) 

i 

c 

m 

o 

■H 

o 

> 

nj 

p 

Q 

0) 

bo 

4-> 

"3 

03 

3 

X> 

r^ 

vO 

X 

CTi 

bO 

r-H 

fn 

a) 

* 

C 

^f 

<D 

X 

<4-l 

03 

O 

2 

X 

<u 

fn 

U 

ol 

ol 

6 

hh 

§ 

U 

3 

3 

UJ 

t/i 

-* 

<0 

1— I 

XI 

ol 

E- 

"3 

C 
03 

c 

nj 
bo 
!h 
o 


3 
X) 

-d 

<D 

>H 

3 
o3 


X 

bO 
U 
0> 

C 
LU 

T3 
CD 


•3 
Cu 


-3 
03 

U 


•3 


H 


o 

O 


O 

o 


o 
o 
o 


o 
o 


o 
o 
o 


o 
o 


o        o 

00  r- 

o       o 


o 

LO 

o 


LO 

o 
o 


CTl 

o 


LO 

o 
o 


o 
o 


LO 

r- 

o 

«* 

CTl 

vO 

O 

o 

LO 

o 

CTl 

o 

CTl 

cti 

O 

r-H 

O 

o 

o 

LO 

o 


o 
o 
o 


to 

CTl 
O 


O 

r-H 

o 


o 


o 
o 


o 

LO 

L0 

tO 

tO 

to 

o 

O 

o 

cti 

CTi 


LO 

o 


o 

00 


00 


o 
o 
o 


o 
o 
o 


o 
o 
o 


o 
to 

CM 


o 
o 


to        o 

o        o 

i— I  (N 


to 

00 


o 

00 


o 


to 


00 

o 


CTi 

o 


r-H 

*tf 

\£) 

LO 

\D 

r» 

\D 

o 

«* 

«* 

<* 

L0 

en 

LO 

to 

oo 

r- 

r-~ 

LO 

00 

t~- 

rl- 

O 

■*D 

■-o 

r-- 

00 

CTl 

CTl 

Cti 

to 

to 

to 

to 

to 

to 

to 

to 

to 

r~ 

to 

«* 

^O 

o 

t"- 

to 

O 

00 

r- 

CTi 

CTl 

00 

CTi 

o 

CTi 

CTl 

^t- 

o 

00 

o 

O 

o 

O 

r-H 

O 

O 

o 

o 

o 

00 
LO 

o 


o 


CTi 
CT) 


LO 


o 

00 


CM 


O 

I— I 

o 

o 


o 
o 
o 


o 
o 


to 

CTi 
O 


o 

o 


O 

o 


o 

1 — I 

o 


o 
o 

o 


oo 
o 


o 
o 
o 


o 
o 
o 


o 

o 

o 

LO 

00 

00 

r-~ 

r-» 

O 

o 

o 

o 

CM  \D 

CTi  00 

O  O 


o 

o 

o 

t~- 

** 

o 

vO 

O 

o 

o 

o 

o 

o 

o 

o 

o 

to 

to 

CTi 

CTi 

vO 

r-~ 

o 

o 

o 

o 

o 

o 

o 

to 

r-~ 

CTi 

CTi 

r-~ 

to 

o 

o 

o 

o 

1— 1 

to 

o 

o 

^o 

md 

00 

LO 

«* 

to 

o 

oo 

r~ 

^D 

LO 

LO 

L0 

LO 

-* 

T* 

*t 

CTl 

*t 

CM 

00 

00 

00 

to 

to 

to 

CTi 
O 


O 

o 
o 


o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

r-H 

to 

LO 

r-~ 

CTi 

(—1 

to 

LO 

r^ 

CTl 

r-H 

to 

r-H 

O 

o 

o 

o 

o 

rsi 

CM 

o 

47 


0.50  I — I — I — •— ^— » — • — *■ 


0.55--. 


0.60 


0.65 --t 


H — I — i — I — I 1 — I 1 1 — I 1 1 1 f- 


0100      0300      0500      0700     0900      1100      1300      1500      1700      1900      2100     2300      0100 


Figure  24.    The  computed  temperature  field  up  to  0.50  bar  for  May  4,  1967  at  Davis,  California. 
First  solid  line  is  280°K  with  a  2°  interval. 
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Figure  25.    The  computed  temperature  field  up  to  2000  meters  for  May  4,  1967  at  Davis,  California. 
First  solid  line  is  290°K  with  a  2°  interval. 
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Figure  26.    Temperoture  profiles  0100  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  27.    Temperature  profiles  0700  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  28.    Temperature  profiles  1300  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  29.    Temperature  profiles  1900  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  30.    Temperature  profiles  0100  PST  on  May  5,  1967  at  Davis/California. 
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Figure  31.    The  computed  potential  temperature  field  up  to  0.50  bar  for  May  4,  1967  at  Davis, 
California. 
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Figure  32.    The  computed  potential  temperature  field  up  to  2000  meters  for  May  4,  1967  at  Davis, 
California. 
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Figure  33.    The  computed  specific  humidity  (gm/kgm)  field  up  to  0.50  bar  for  May  4,  1967  at  Davis, 
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Figure  34.    The  computed  specific  humidity  (gm/kgm)  field  up  to  2000  meters  for  May  4,  1967  at 
Davis,  California. 
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Figure  35.    Specific  humidity  profiles  0100  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  36.    Specific  humidity  profiles  0700  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  37.    Specific  humidity  profiles  1300  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  38.    Specific  humidity  profiles  1900  PST  on  May  4,  1967  at  Davis,  California. 
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Figure  39.    Specific  humidity  profiles  0100  PST  on  May  5,  1967  at  Davis,  California. 
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Figure  40.    The  computed  relative  humidity  field  up  to  0.50  bar  for  May  4,  1967  at  Davis,  California. 
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Figure  41.    The  computed  relative  humidity  field  up  to  2000  meters  for  May  A,  1967  at  Davis, 
California. 
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4.4  Calculations  Using  Different  Boundary  Conditions 

In  the  O'Neill  calculation  the  specific  humidity  at  the  earth's  surface  was  specified 
as  a  function  of  time,  and  in  the  Davis  calculation  the  latent  heat  flux  at  the  earth's 
surface  was  specified  as  a  function  of  time.   Calculations  were  also  made  where  (1)  the  rel- 
ative humidity  at  the  earth's  surface  was  held  constant  and  (2)  where  the  specific  humidity 
at  the  earth's  surface  was  held  constant. 

Using  a  constant  relative  humidity  of  80  percent  at  the  earth's  surface  in  the  Davis 
calculation  resulted  in  predicted  temperatures  which  were  slightly  less  than  those  resulting 
from  the  calculation  described  above.   The  lower  temperatures  are  attributed  to  a  larger 
latent  heat  flux  at  the  earth's  surface.  The  maximum  increase  in  the  predicted  specific 
humidity  was  1.0  gm/kgm,  and  the  maximum  decrease  in  temperature  was  1.5°C.   The  slight  var- 
iations between  the  two  calculations  are  due  primarily  to  the  fact  that  the  diurnal  tempera- 
ture variation  was  only  about  5°C  near  the  earth's  surface,  and  the  observed  relative  humid- 
ity did  not  vary  greatly. 

The  results  were  not  as  satisfactory  when  the  relative  humidity  was  held  constant  in  the 
O'Neill  calculation.  The  diurnal  temperature  variations  resulted  in  variations  in  the  spe- 
cific humidity  which  were  out  of  phase  with  the  observed  variations.   The  predicted  specific 
humidity  increased  in  the  afternoon  and  decreased  between  midnight  and  dawn  as  the  tempera- 
ture changed,  whereas  the  observed  specific  humidity  decreased  in  the  afternoon  and  increased 
between  midnight  and  dawn. 

Variations  in  the  specific  humidity  such  as  those  observed  at  O'Neill  are  typical  of 
the  effect  of  stomatal  closure.   Stomatal  closure  reduces  the  amount  of  water  vapor  released 
from  the  plants  by  transpiration,  resulting  in  a  lower  specific  humidity  in  the  afternoon. 
This  demonstrates  the  need  to  include  plant  physiological  processes  in  a  moisture  balance 
calculation  at  the  earth's  surface. 

When  a  constant  specific  humidity  of  6.70  gm/kgm  at  the  earth's  surface  was  used  in  the 
Davis  calculation,  the  predicted  temperatures  were  slightly  greater  than  those  resulting 
from  specifying  the  latent  heat  flux.  The  predicted  latent  heat  flux  generally  ranged  from 
50  to  90  percent  of  the  values  for  the  other  calculation,  which  accounted  for  the  higher 
predicted  temperatures.   The  maximum  temperature  difference  was  1.5°C,  and  the  maximum  dif- 
ference in  the  specific  humidities  was  1.3  gm/kgm.   The  values  of  the  predicted  specific 
humidity  were  smaller  for  this  calculation. 

When  a  constant  specific  humidity  of  12.7  gm/kgm  was  used  in  the  O'Neill  calculation, 
the  predicted  temperatures  varied  from  1.7°C  less  to  2.0  degrees  more  than  those  resulting 
from  specifying  the  specific  humidity  at  the  earth's  surface.   At  the  end  of  the  calculation 
the  predicted  temperatures  were  smaller  by  about  0.2°C.  The  latent  heat  flux  was  greater 
for  this  calculation  during  the  first  part  of  the  calculation,  and  it  was  smaller  during  the 
last  part  of  the  calculation.  The  predicted  moisture  field  showed  poor  agreement  with  the 
observed  moisture  field  because  the  large  variations  in  the  observed  moisture  field  could 
not  be  predicted  by  holding  the  specific  humidity  constant  at  the  earth's  surface. 

4. 5  Advection 

Estoque  attributed  some  of  the  discrepancy  between  his  calculated  temperature  field  and 
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the  observed  temperature  field  to  the  effects  of  advection,  but  he  did  not  estimate  the  mag- 
nitude of  the  horizontal  temperature  gradient  necessary  to  account  for  the  discrepancy. 
The  O'Neill  site  was  chosen  because  it  was  thought  to  be  nearly  ideal.   The  surrounding  ter- 
rain is  flat  with  a  uniform  ground  cover,  which  should  minimize  the  effects  of  advection, 
yet  the  magnitude  of  the  advection  has  not  been  determined. 

The  magnitude  of  the  horizontal  temperature  gradients  required  to  bring  the  predicted 
temperature  profile  into  agreement  with  the  observed  temperature  profile  were  calculated  for 
the  six-hour  period  between  1230  and  1830  CST  on  August  24.   The  horizontal  temperature  gra- 
dients were  calculated  as  a  function  of  height  using  a  trial  and  error  solution.   The  hor- 
izontal temperature  gradients  for  the  region  below  200  meters  are  summarized  in  Table  5.   The 
x  direction  is  defined  as  the  direction  of  the  wind  vector. 

Table  5:   The  horizontal  temperature  gradient  at  100  meters  on 
August  24,  1953  at  O'Neill,  Nebraska 

Time  LST  |I  -u|I 

3x  3x 


1230  -  1430         0.000°C/km  0.0  °C/hr 

1430  -  1630         0.008°C/km         -0.43°C/hr 
1630  -  1830         0.019°C/km         -0.88°C/hr 

A  rough  estimate  for  the  actual  horizontal  temperature  gradient  may  be  made  from  the 
temperature  data  for  the  region  surrounding  the  test  site,  which  is  shown  in  Figure  B.5.3.a 
of  the  O'Neill  data  (Lettau  and  Davidson,  1953).   According  to  these  data,  the  horizontal 
temperature  gradient  near  the  earth's  surface  at  1530  CST  on  August  24  could  range  from  about 
0.005°C/km  to  0.010°C/km  with  cooler  air  upwind  from  the  site.   Since  the  calculated  temper- 
ature gradient  lies  within  this  range  of  values,  this  tends  to  support  the  argument  that  the 
discrepancies  between  the  observed  and  predicted  temperatures  can  be  attributed  to  advection 
effects.   This  also  tends  to  verify  the  accuracy  of  the  other  processes  in  the  calculation 
because  it  was  not  necessary  to  account  for  the  discrepancies  by  an  unreasonably  large  hori- 
zontal temperature  gradient. 

Although  the  horizontal  temperature  gradients  were  small,  they  were  sufficient  to  cause 
a  temperature  change  of  more  than  2.5°C  over  a  four-hour  period  because  of  the  strong  wind. 
This  is  a  significant  effect,  which  cannot  be  neglected. 

The  horizontal  temperature  gradients  were  also  estimated  for  the  Davis  calculation  over 
a  similar  six-hour  period.   Beginning  with  the  observed  temperature  and  moisture  profiles  at 
1300  PST  on  May  4,  1967,  the  horizontal  temperature  gradients  required  to  maintain  agreement 
between  the  observed  and  predicted  temperature  gradients  were  calculated.   The  results  are 
summarized  in  Table  6. 
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Table  6:   The  horizontal  temperature  gradient  at  100  meters  on 
May  4,  1967  at  Davis  California 


Time  LST 

3T 
9x 

-uJl 

3x 

1300  -  1430 

0.00°C/km 

0.00°C/hr 

1430  -  1530 

0.01°C/km 

-0.20°C/hr 

1530  -  1630 

0.09°C/km 

-1.40°C/hr 

1630  -  1730 

0.05°C/km 

-0.99°C/hr 

1730  -  1900 

0.03°C/km 

-0.81°C/hr 

Fosberg  and  Schroeder  (1966)  analyzed  marine  air  penetration  in  central  California  dur- 
ing the  summers  of  1961  and  1962.   They  observed  a  temperature  difference  of  5.8°C  between 
Davis  and  a  point  in  the  Sacramento  Delta  about  40  miles  upwind  of  Davis  at  1400  PST  on 
August  15,  1962.   The  temperatures  were  cooler  in  the  Delta,  which  lies  to  the  south-west 
of  Davis.  The  average  horizontal  temperature  gradient  was  about  0.09°C/km,  which  agrees 
well  with  the  values  in  Table  6,  bearing  in  mind  that  the  data  do  not  correspond  to  the 
same  day. 

The  horizontal  temperature  gradients  presented  in  Table  6  accounted  for  a  temperature 
change  of  about  4°C  over  a  period  of  4-1/2  hours.  Although  the  temperature  gradients  were 
larger  for  the  Davis  calculation  than  for  the  O'Neill  calculation,  the  resulting  temperature 
changes  were  not  significantly  different.   This  demonstrates  the  importance  of  the  combina- 
tion of  wind  speed  and  horizontal  temperature  gradient  when  considering  the  effect  of 
advection. 

The  rate  of  change  of  temperature  due  to  advection  is  less  near  the  earth's  surface  where 
the  winds  are  lighter  than  it  is  at  higher  elevations.  However, significant  changes  may  occur 
in  the  energy  transfer  processes  near  the  earth's  surface  as  a  result  of  advection  at  higher 
levels.  This  is  shown  in  Tables  7  and  8.  The  rates  of  change  of  temperature  due  to  various 
energy  transfer  processes,  not  including  advection,  are  shown  in  Table  7  for  the  Davis  cal- 
culation at  1700  on  May  4,  1967.   A  similar  set  of  data  is  shown  in  Table  8  corresponding  to 
the  same  time,  but  in  this  case  advection  is  included.  The  computed  horizontal  temperature 
gradient  equaled  0.05°C/km  up  to  a  height  of  600  meters,  decreased  to  zero  at  1000  meters, 
and  was  zero  above  this  height. 

Although  the  rate  of  change  of  temperature  due  to  advection  varies  a  great  deal  in  this 
region,  the  total  rate  of  cooling  increased  nearly  uniformly  when  advection  was  included. 
This  can  be  explained  by  looking  at  the  rate  of  change  of  temperature  due  to  convection, 
which  is  where  the  most  significant  changes  occur.  According  to  Table  7  where  advection  is 
not  included,  heat  is  moved  downward  by  convection,  thus  warming  the  region  near  the  earth's 
surface  and  cooling  the  region  at  higher  levels.  When  advection  is  included,  as  shown  in 
Table  8,  heat  is  moved  upward  into  the  area  being  cooled  most  rapidly  by  advection.   The 
result  is  that  the  region  near  the  earth's  surface  is  cooled  by  convection  instead  of  being 
warmed  by  convection  as  was  the  case  before. 


59 


c 

o 

•H 

■)-> 

o 

ID 

> 

TJ 

< 

b£> 

c 

•H 

T) 

P 

i— 1 

u 

X 

d 

LU 

■  H 

C 

(D 

'H 

u 

o 

0) 

4-i 

X 

■H 

&. 

i— 1 

10 

rt 

o 

U 

e 

+j 

* 

< 

i/i 

•H 

h 

> 

<L> 

Cfl 

2 

Q 

O 

J 

+-» 

cti 

<U 

x: 

t-^ 

■M 

s0 

o> 

0) 

M 

X 

TJ 

a} 

3 

s 

CO 

c 

X 

o 

W) 

u 

H 

CD 

CO 

c 

D. 

UJ 

O 

0) 

o 

x: 

t-^ 

H 

t— 1 

r». 

<u 

^H 

x> 

tti 

H 

(-1 

o  ~^ 

m  u 

E-  o 


Q  -C 
E-  o 


>    M 

Z  ,C 

O  ■>«. 
u  u 

E-  o 

Q  »— 


.J  M 
<  X 
E-i  --* 
O  U 
H  o 


Ho 


m  e 


Ol 

r— I 
O 

oo 
o 

1— 1 

o 

f— 1 

o 

.— 1 

o 

o 

i-H 

o 

t— 1 

o 

i— I 

o 

o 

i— 1 

o 

i— 1 
O 

i— 1 
O 

r-H 
O 

o 

i — i 

o 

i-H 

o 

1— 1 

o 

00 

r-H 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o 

o 

O 

O 

o 

o 

o 

o 

o 

CN       tO 

o     o 


r-~     to 

vO       00 

o     o 


© 


t~-       00 

00       00 

cn     cn 


CTl 

o 


1^ 


LO       00       CTl       rf 
vD       vD       v£)       1^ 

o     o     o     o 


CN  "3- 

oo     en 
o     o 


1— 1 

LO 

\Q 

CN) 

o 

to 

LO 

o 

CTl 

vO 

LO 

vO 

CN 

to 

tO 

LO 

00 

^i- 

en 

vO 

(XI 

to 

CTl 

o 

r*. 

CT) 

CM 

-T 

LO 

^D 

vO 

LO 

LO 

to 

i — ( 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

CTl 

00 

CN 

to 

•■o 

vO 

o 

h- 

CTl 

a, 

o 

t— 1 

t— 1 

CN 

■*t 

o 


to 


to 


o 


CTl 
O 


CN 

O 


O 


00 
CN 


o 


O 
00 


O 
00 


Ol 


oo 


CN 


en 

00 
CN 


o 

CT, 
CN 


o 

CT. 
CN 


o 

CTl 
CN 


CTl 
CN 


CTl 
CN 


CTl 
CN 


to 


to 

o 

Ol 

CN 
LO 

CN 

to 

i-H 

r-H 
LO 

O 

00 

LO 

t— t 

CN 

LO 

CD 
i— I 

to 

to 

00 

o 

LO 

00 

LO 

CTl 

to 

O 

to 

(N 
<N 

t~- 

LO 

to 

CN 

to 

LO 

CT>i— (LOOvO^CNCNCNi— It— It— lOOCTlOOvDtOr~» 
tOtOCNCNt-Hi— It-Ht— It— IrHt— It-H,— It-IOOOOCTI 

■-Ir-lT-lr-li-ll-ll-lr-ll-tl-il-il-ll-ll-llHl-t.f-Hl-lO 

CTlCTlOlCTlOlOlOlOlCTlOlOlOlOlOlCTlOlOlOlOl 
CNCNCNCNCNCNCNCNCNCNCNCNCNCNCNCNCNCNCN 


00 

CTl 


Ol       Ol      CTl       CTl       CTl       CTl 


o 

LO 


LO 

CTl 


Ol 
CN 


o 

o 


60 


o 

•H 

*-> 

o 

a> 

> 

Ti 

< 

00 

c 

•H 

TJ 

3 

■— i 

(J 

C 

ci 

HH 

•H 

c 

<D 

R 

h 

o 

<D 

4h 

X 

•H 

p< 

. — l 

w 

Rj 

o 

U 

6 

+j 

* 

< 

t/> 

■H 

Sh 

> 

(U 

nt 

3 

Q 

O 

-J 

*J 

cd 

0) 

x; 

r- 

♦j 

vD 

a, 

a> 

W> 

>, 

T) 

crt 

3 

S 

CO 

P! 

>, 

O 

M 

U 

H 

a> 

CO 

c 

Cu 

OJ 

O 

<D 

o 

X 

r-^ 

H 

(— 1 

00 

<u 

t— 1 

£> 

crt 

H 

o  --. 
CO  u 
H  o 


a  x: 
2  u 

H  o 


>    X 

Q  --. 
<    U 
H  o 
Q   "— 


O  --s. 

U  U 

H  o 

Q  ^ 


<  X 

H  --. 

O  U 

H  o 


CD   O 


H  o 


M    6 


to 

O 

CM 
CM 
O 

i— i 
O 

o 

CM 

o 

at 
r— i 
o 

at 

i—i 
o 

at 
i— i 
o 

oo 

i— ( 
o 

00 

1—1 

o 

00 
1— 1 

o 

00 

i— 1 

o 

oo 

1— I 

o 

00 

1— 1 

o 

i-H 

o 

i-H 

o 

00 

t— l 

o 

00 
O 

oo 
o 

at 
i—i 
o 

O 

O 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

00 


at 

LO 

o 


cr. 
co 


LO         CM         LO         LO 

LO       \£)       v£>       vO 
O       O       O       O 


o 
to 


CTi 


CO 
CO 


CTI 
O 


-3- 

co 
o 


en 
o 


LO 


00 
O 


CO 

O 

CO 

at 

to 

LO 

CO 

o 

to 

\Q 

to 

OI 

to 

00 

o 


CTl 

CT. 


cr.      oo      h~ 


a-, 

i-H 

00 

LO 

to 

CM 

CM 

CM 

o 

o 

at 

o 

CM 

LO 

CTi 

<* 

vD 

LO 

to 

to 

CM 

i— I 

O 

O 

vo     r~»     oo     oo 
•^-     r^     to     rt- 

^t      Tt       Tt      to 


IO 

at 

CM 

o 

o 

vO 

i-H 

i-H 

\0 

\0 

r-. 

O 

00 

l-- 

\o 

vO 

L0 

to 

i — i 

o 

O 

o 

r— 1 

CM 

to 

«tf 

LO 

CT. 

to 

CO 

tt 

to 

-t 

vO 

vD 

LO 

<* 

LO 

*fr 

00 

CO 

CT. 

CM 

"■* 

LO 

vD 

CO 

CTi 

at 

at 

CTi 

<* 

o 

LO 

CT. 

to 

CO 

CO 

CTi 

CTl 

CT. 

O 

o 

a. 

CTi 

CT. 

CT. 

o 

o 

to 
o 


CTi 
O 


O        i-H        i-H 


CT. 
CO 


LO 

co 


o 

00 


o 

CO 


o 
oo 


o 
co 


00        00        00        00        00        00 


co 


00 


CO 

co 


00 
00 


co 

CO 
CM 


00 

CO 
CM 


00 
00 


co 
co 

CM 


OO 
00 
CM 


00 

00 


co 

CO 


00  00 
00  00 
CM        CM 


o 

CTl 


CO 

co 


o 
to 


CTi 


co 

co 

CM 


CM 

oo 

i-H 

o 

LO 

LO 

i— I 

o 

r— 1 

LO 

CTl 

LO 

i-H 

CM 

i-H 

i— 1 
LO 

00 

LO 

L0 

00 

<* 

00 

to 

CTi 

>* 

co 

CM 

LO 
00 
CM 

00 
CM 

00 
CM 

00 
CM 

00 
00 
CM 

00 
00 
CM 

oo 

00 
CM 

00 
00 
CM 

CTi 
00 
CM 

CTi 
00 
CM 

at 

00 
CM 

CTl 
00 
CM 

CTl 
00 
CM 

CTi 
00 
CM 

CTi 
00 
CM 

at 

00 
CM 

CTi 
00 
CM 

CTl 
00 
CM 

i-H 

CTi 

r-» 

CM 

t~- 

O 

CM 

CM 

to 

1— 1 

o 

00 

LO 

CM 

i— 1 

to 

to 

o 

o 

i — i 

O 

CM 

CM 

to 

en 

LO 

Tf- 

I—i 

LO 

vO 

vO 

i-H 

LO 

•* 

LO 

vi) 

LO 

o 

f» 

vO 

r» 

vO 

CTi 

to 

t~» 

00 

LO 

\0 

i-H 

CTi 

O 

CM 

^0 

,_, 

t^ 

rf 

CM 

-t- 

to 

-* 

r^ 

r— 1 

r^ 

to 

o 

00 

vO 

LO 

to 

to 

CM 

i-H 

1—1 

m 

■5* 

to 

CM 

CM 

i-H 

i—l 

i-H 

61 


5.0  CONCLUSIONS 

One  of  the  primary  goals  of  this  work  was  the  development  of  a  numerical  model  which  is 

more  versatile  and  as  accurate  or  more  accurate  than  existing  numerical  models  in  predicting 

temperature  and  moisture  fields  in  the  lower  atmosphere.   The  O'Neill  calculation  demonstrated 

that  the  present  numerical  model  is  as  accurate  or  more  accurate  than  Estoque's  model.   This 

is  significant  because  the  present  model  was  designed  to  apply  to  situations  with  a  broad 

i 
range  of  weather  conditions,  whereas  Estoque's  model  was  designed  primarily  to  handle  the 

special  conditions  present  at  O'Neill. 

The  Davis  calculation  demonstrated  the  ability  of  the  numerical  model  to  predict  tem- 
perature and  humidity  fields  under  conditions  of  light  winds.   Accurate  wind  data  must  be 
available  for  the  calculation  of  accurate  sensible  heat  fluxes  near  the  earth's  surface. 
This  seemed  to  be  particularly  important  under  conditions  of  light  wind  and  strong  thermo- 
dynamic stability  near  the  earth's  surface  which  make  the  calculation  of  accurate  convective 
heat  fluxes  more  difficult. 

Comparison  of  calculated  values  of  the  constituents  of  the  energy  budget  at  the  earth's 
surface  with  measured  values  for  both  the  Davis  and  O'Neill  calculations  showed  the  terres- 
trial and  solar  radiation  calculations  to  be  very  accurate.   Satisfactory  results  were  ob- 
tained for  the  water  vapor  calculation  by  either  specifying  the  latent  heat  flux  or  the 
specific  humidity  at  the  earth's  surface  as  functions  of  time.   Improvements  must  be  made 
in  the  water  vapor  calculation,  however,  if  the  formation  of  fog  or  stratus  is  to  be  pre- 
dicted accurately. 

Advection  is  a  very  significant  energy  transfer  process  which  should  not  be  neglected 
in  future  numerical  models.   In  the  past  the  assumption  of  horizontal  homogeneity  was  made 
primarily  to  allow  the  use  of  a  one  dimensional  model,  which  requires  less  calculating  time 
and  less  computer  storage  space.   With  the  advent  of  faster  computers  with  greater  storage, 
the  assumption  of  horizontal  homogeneity  should  not  be  necessary.  A  multi-dimensional  model 
should  lead  to  more  accurate  calculation  of  the  effects  of  advection,  which  would  be  an  im- 
portant step  toward  further  improvement  of  weather  prediction  by  numerical  simulation.  One- 
dimensional  models  such  as  this  one  may  prove  useful  for  qualitative  studies  or  special  case 
studies,  and  their  importance  should  not  be  overlooked. 

Advection  has  a  significant  effect  on  the  convective  heat  transfer  near  the  earth's 
surface  even  though  the  rate  of  change  of  temperature  due  to  advection  may  be  small  in  this 
region.   Advection  should  not  be  neglected  in  specialized  numerical  models  which  model  only 
a  few  meters  of  the  atmosphere,  such  as  might  be  used  for  the  prediction  of  radiation  fog, 
because  of  its  effect  on  the  other  energy  transfer  processes. 

The  models  for  the  eddy  diffusion  coefficient  suggested  by  Deardorff  (1967)  and  Craw- 
ford (1965)  and  the  terrestrial  radiation  calculation  suggested  by  Crawford  (1965)  were  used 
successfully  in  this  model  in  modified  form.  Hopefully,  the  results  presented  here  will  lead 
to  further  improvements  in  these  models. 

The  model  for  the  eddy  diffusion  coefficient  could  be  improved  for  conditions  of  strong 
thermodynamic  stability.   This  may  be  done  most  easily  by  modifying  <t>(R.)  for  Richardson 
numbers  greater  than  0.1.   The  model  for  the  eddy  diffusion  coefficient  may  also  be  improved 
by  allowing  z  and  z~   to  vary  during  the  calculation. 
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The  longwave  and  solar  radiation  calculations  can  be  improved  by  including  fog  or 
stratus  in  the  calculation.   Eventually  precipitation  and  the  effect  of  evaporation  of  the 
falling  droplets  should  be  included  in  the  model. 

The  moisture  calculation  could  be  improved  by  using  a  moisture  balance  calculation  at 
the  earth's  surface.   Plant  physiological  processes  should  also  be  included  in  the  calcul- 
ation.  These  improvements  would  eliminate  the  need  to  specify  any  of  the  moisture  con- 
ditions at  the  earth's  surface. 

The  soil  heat  flux  calculation  could  be  improved  by  allowing  the  soil  density  and  thermal 
conductivity  to  vary  with  depth  and  time.   Different  composition  of  the  soil  and  variations 
in  the  water  content  would  cause  the  values  of  these  parameters  to  vary. 

Mating  this  model  to  a  general  circulation  model  may  lead  to  improvements  in  both 
models.   The  general  circulation  calculation  would  be  improved  by  including  detailed  weather 
conditions  near  the  earth's  surface  which  would  not  otherwise  be  apparent  in  the  calculation. 
The  present  model  would  be  improved  because  conditions  in  the  upper  atmosphere  would  be  al- 
lowed to  influence  conditions  in  the  lower  atmosphere.   The  conditions  at  the  upper  boundary 
of  the  atmospheric  region  of  the  present  model  would  be  modified  to  allow  fluxes  of  heat 
and  water  vapor  across  the  boundary. 

It  is  hoped  that  this  work  will  lead  to  further  attempts  at  forecasting  the  weather  by 
numerical  simulation.  Only  by  such  attempts  will  accurate  prediction  of  the  weather  become 
a  reality. 
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